Nothing
## ----echo=FALSE, results='hide'-----------------------------------------------
# This should be set to FALSE in case you want to redo the sampling as well.
# Default is TRUE to save CRAN some time.
usePreCalcResults <- TRUE
# Default here is FALSE because the object created is rather large
usePreCalcResultsExtra <- FALSE
## ----setup, include=FALSE, cache=FALSE------------------------------
knitr::render_sweave()
knitr::opts_chunk$set(prompt = TRUE,
fig.show = "hide",
fig.keep = "high",
warning = FALSE,
error = FALSE,
dev = c("png", "pdf"),
dpi = 90,
message = FALSE,
echo = FALSE,
eval=TRUE,
#cache = TRUE,
fig.path = "Figures/jss2014-",
tidy = FALSE)
base::options(continue = "+ ", prompt = "R> ", width = 70,
useFancyQuotes = FALSE)
plotevery <- 1L # can be used for trace plots to make vignette smaller
## ----usd1, echo=TRUE, fig.width=11, fig.height=6--------------------
set.seed(123)
library("stochvol")
data("exrates")
ret <- logret(exrates$USD, demean = TRUE)
par(mfrow = c(2, 1), mar = c(1.9, 1.9, 1.9, 0.5), mgp = c(2, 0.6, 0))
plot(exrates$date, exrates$USD, type = "l",
main = "Price of 1 EUR in USD")
plot(exrates$date[-1], ret, type = "l", main = "Demeaned log returns")
## ----usd2, echo=TRUE, fig.width=11, fig.height=6--------------------
sim <- svsim(500, mu = -9, phi = 0.99, sigma = 0.1)
par(mfrow = c(2, 1))
plot(sim)
## ----svsample, eval=TRUE, echo=TRUE---------------------------------
ret <- ret[1:200]
res <- svsample(ret, priormu = c(-10, 1), priorphi = c(20, 1.1),
priorsigma = 0.1, thin = 10)
## ----summary, echo=TRUE---------------------------------------------
summary(res, showlatent = FALSE)
## ----paropar, include=FALSE-----------------------------------------
par(mfrow = c(1, 1))
## ----usd3, echo=TRUE, fig.width=11, fig.height=7--------------------
volplot(res, forecast = 100, dates = exrates$date[seq_along(ret)])
## ----usd4, echo=TRUE, fig.width=11, fig.height=7--------------------
res <- updatesummary(res, quantiles = c(0.01, 0.1, 0.5, 0.9, 0.99))
volplot(res, forecast = 100, dates = exrates$date[seq_along(ret)])
## ----usd5, echo=TRUE, fig.width=8, fig.height=5---------------------
par(mfrow = c(3, 1))
paratraceplot(res)
## ----usd6, echo=TRUE, fig.width=8, fig.height=3.5-------------------
par(mfrow = c(1, 3))
paradensplot(res, showobs = FALSE)
## ----usd7, echo=TRUE, fig.width=9, fig.height=7---------------------
plot(res, showobs = FALSE)
## ----usd8, echo=TRUE, fig.width=9, fig.height=7---------------------
myresid <- resid(res)
plot(myresid, ret)
## ----beta1, echo=TRUE-----------------------------------------------
set.seed(123456)
n <- 200
beta.true <- c(0.1, 0.5)
sigma.true <- 0.01
X <- matrix(c(rep(1, n), rnorm(n, sd = sigma.true)), nrow = n)
y <- rnorm(n, X %*% beta.true, sigma.true)
## ----betaprior, echo=TRUE-------------------------------------------
burnin <- 100
draws <- 1000
b0 <- matrix(c(0, 0), nrow = ncol(X))
B0inv <- diag(c(10^-10, 10^-10))
c0 <- 0.001
C0 <- 0.001
## ----betaprecalc, echo=TRUE-----------------------------------------
p <- ncol(X)
preCov <- solve(crossprod(X) + B0inv)
preMean <- preCov %*% (crossprod(X, y) + B0inv %*% b0)
preDf <- c0 + n/2 + p/2
## ----echo=TRUE------------------------------------------------------
draws1 <- matrix(NA_real_, nrow = draws, ncol = p + 1)
colnames(draws1) <- c(paste("beta", 0:(p-1), sep = "_"), "sigma")
sigma2draw <- 1
## ----betaloop, echo=TRUE--------------------------------------------
for (i in -(burnin-1):draws) {
betadraw <- as.numeric(mvtnorm::rmvnorm(1, preMean,
sigma2draw * preCov))
tmp <- C0 + 0.5 * (crossprod(y - X %*% betadraw) +
crossprod((betadraw - b0), B0inv) %*% (betadraw - b0))
sigma2draw <- 1 / rgamma(1, preDf, rate = tmp)
if (i > 0) draws1[i, ] <- c(betadraw, sqrt(sigma2draw))
}
## ----colmeansbeta, echo=TRUE----------------------------------------
colMeans(draws1)
## ----draws1, eval=TRUE, echo=TRUE, fig.width=7.5, fig.height=6------
plot(coda::mcmc(draws1))
## ----echo=TRUE------------------------------------------------------
mu.true <- log(sigma.true^2)
phi.true <- 0.97
vv.true <- 0.3
simresid <- svsim(n, mu = mu.true, phi = phi.true, sigma = vv.true)
y <- X %*% beta.true + simresid$y
## ----echo=TRUE------------------------------------------------------
draws <- 5000
burnin <- 500
thinning <- 10
priors <- specify_priors(
mu = sv_normal(-10, 2),
phi = sv_beta(20, 1.5),
sigma2 = sv_gamma(0.5, 0.5)
)
## ----echo=TRUE------------------------------------------------------
draws2 <- matrix(NA_real_, nrow = floor(draws / thinning),
ncol = 3 + n + p)
colnames(draws2) <- c("mu", "phi", "sigma",
paste("beta", 0:(p-1), sep = "_"), paste("h", 1:n, sep = "_"))
betadraw <- c(0, 0)
paradraw <- list(mu = -10, phi = 0.9, sigma = 0.2)
latentdraw <- rep(-10, n)
paranames <- names(paradraw)
## ----eval=FALSE, echo=TRUE------------------------------------------
# for (i in -(burnin-1):draws) {
# ytilde <- y - X %*% betadraw
# svdraw <- svsample_fast_cpp(ytilde, startpara = paradraw,
# startlatent = latentdraw, priorspec = priors)
# paradraw <- svdraw$para
# latentdraw <- drop(svdraw$latent)
# normalizer <- as.numeric(exp(-latentdraw / 2))
# Xnew <- X * normalizer
# ynew <- y * normalizer
# Sigma <- solve(crossprod(Xnew) + B0inv)
# mu <- Sigma %*% (crossprod(Xnew, ynew) + B0inv %*% b0)
# betadraw <- as.numeric(mvtnorm::rmvnorm(1, mu, Sigma))
# if (i > 0 && i %% thinning == 0) {
# draws2[i/thinning, 1:3] <- drop(paradraw)[paranames]
# draws2[i/thinning, 4:5] <- betadraw
# draws2[i/thinning, 6:(n+5)] <- latentdraw
# }
# }
## ----longerbeta-----------------------------------------------------
if (usePreCalcResults) {
load("vignette_sampling_draws2.RData")
} else {
for (i in -(burnin-1):draws) {
# draw latent volatilities and AR-parameters:
ytilde <- y - X %*% betadraw
svdraw <- svsample_fast_cpp(ytilde, startpara = paradraw,
startlatent = latentdraw, priorspec = priors)
paradraw <- svdraw$para
latentdraw <- drop(svdraw$latent)
# draw the betas:
normalizer <- as.numeric(exp(-latentdraw/2))
Xnew <- X * normalizer
ynew <- y * normalizer
Sigma <- solve(crossprod(Xnew) + B0inv)
mu <- Sigma %*% (crossprod(Xnew, ynew) + B0inv %*% b0)
betadraw <- as.numeric(mvtnorm::rmvnorm(1, mu, Sigma))
# store the results:
if (i > 0 & i %% thinning == 0) {
draws2[i/thinning, 1:3] <- drop(paradraw)[paranames]
draws2[i/thinning, 4:5] <- betadraw
draws2[i/thinning, 6:(n+5)] <- latentdraw
}
}
draws2selection <- draws2[,4:8]
save(draws2selection, file = 'vignette_sampling_draws2.RData')
}
## ----betaplot2, echo=TRUE, eval=FALSE-------------------------------
# plot(coda::mcmc(draws2[, 4:8]))
## ----hetero, fig.width=7.5, fig.height=8----------------------------
par(mar = c(3.1, 1.8, 1.9, .5), mgp = c(1.8, .6, 0))
plot(coda::mcmc(draws2selection[seq(1L, nrow(draws2selection), by = plotevery), 1:4]))
## ----colmeans2, eval=FALSE, echo=TRUE-------------------------------
# colMeans(draws2[, 4:8])
## ----echo=FALSE, eval=TRUE------------------------------------------
colMeans(draws2selection)
## ----scatter, echo=TRUE, fig.width=10.2-----------------------------
x <- log(exrates$USD[-length(exrates$USD)])
y <- log(exrates$USD[-1])
X <- matrix(c(rep(1, length(x)), x), nrow = length(x))
par(mfrow=c(1,1), mar = c(2.9, 2.9, 2.7, .5), mgp = c(1.8,.6,0), tcl = -.4)
plot(x,y, xlab=expression(log(p[t])), ylab=expression(log(p[t+1])),
main="Scatterplot of lagged daily log prices of 1 EUR in USD",
col="#00000022", pch=16, cex=1)
abline(0,1)
## ----results='hide'-------------------------------------------------
if (usePreCalcResults) {
load("vignette_sampling_realworld.RData")
} else {
set.seed(34567890)
#configuration parameters
draws <- 100000
thinning <- 100
burnin <- 10000
priormu <- c(-10, 1)
priorphi <- c(20, 1.5)
priorsigma <- .1
p <- 2
n <- length(x)
## design matrix:
X <- matrix(c(rep(1, n), x), nrow=n)
## prior for beta (or beta|sigma in homoskedastic regression):
b0 <- matrix(c(0, 0), nrow=ncol(X))
B0 <- diag(c(10^10, 10^10))
## prior for sigma^2 (only used in homoskedastic regression)
c0 <- 0.001
C0 <- 0.001
B0inv <- solve(B0)
# specify priors in an object
priors <- specify_priors(
mu = sv_normal(-10, 2),
phi = sv_beta(20, 1.5),
sigma2 = sv_gamma(0.5, 0.5)
)
## initialize some space to hold results:
realres <- vector("list", 3)
## AR(1)-SV:
realres[[1]] <- matrix(NA_real_, nrow = floor(draws/thinning), ncol = 3 + n + p)
colnames(realres[[1]]) <- c("mu", "phi", "sigma", paste("beta", 0:(p-1), sep="_"),
paste("h", 1:n, sep="_"))
## some indicators:
paras <- 1:3
betas <- 3+(1:p)
latents <- 3+p+(1:n)
## starting values:
betadraw <- rep(.1, p)
paradraw <- list(mu = -10, phi = 0.8, sigma = 0.2)
latentdraw <- rep(-10, n)
paranames <- names(paradraw)
## sampler:
tim <- system.time(
for (i in -(burnin-1):draws) {
if (i%%1000 == 0) cat("Iteration", i, "done.\n")
# draw latent volatilities and SV-parameters:
ytilde <- y - X %*% betadraw
svdraw <- svsample_fast_cpp(ytilde, startpara = paradraw,
startlatent = latentdraw, priorspec = priors)
paradraw <- svdraw$para
latentdraw <- drop(svdraw$latent)
# draw the betas:
normalizer <- as.numeric(exp(-latentdraw/2))
Xnew <- X * normalizer
ynew <- y * normalizer
Sigma <- solve(crossprod(Xnew) + B0inv)
mu <- Sigma %*% (crossprod(Xnew, ynew) + B0inv %*% b0)
betadraw <- as.numeric(mvtnorm::rmvnorm(1, mu, Sigma))
# store the results:
if (i > 0 & i %% thinning == 0) {
realres[[1]][i/thinning, paras] <- drop(paradraw)[paranames]
realres[[1]][i/thinning, latents] <- latentdraw
realres[[1]][i/thinning, betas] <- betadraw
}
}
)[["elapsed"]]
## AR(1) homoskedastic:
## pre-calculate some values:
preCov <- solve(crossprod(X) + B0inv)
preMean <- preCov %*% (crossprod(X, y) + B0inv %*% b0)
preDf <- c0 + n/2 + p/2
## allocate space:
realres[[2]] <- matrix(as.numeric(NA), nrow=floor(draws/thinning), ncol=p+1)
colnames(realres[[2]]) <- c(paste("beta", 0:(p-1), sep="_"), "sigma")
## starting values:
betadraw <- rep(0, p)
sigma2draw <- var(y)
## sampler:
tim2 <- system.time(
for (i in -(burnin-1):draws) {
if (i%%1000 == 0) cat("Iteration", i, "done.\n")
# draw beta:
betadraw <- as.numeric(mvtnorm::rmvnorm(1, preMean, sigma2draw*preCov))
# draw sigma^2:
tmp <- C0 + .5*(crossprod(y - X%*%betadraw) +
crossprod((betadraw-b0), B0inv) %*% (betadraw-b0))
sigma2draw <- 1/rgamma(1, preDf, rate=tmp)
# store the results:
if (i > 0 & i %% thinning == 0) {
realres[[2]][i/thinning, 1:p] <- betadraw
realres[[2]][i/thinning, p+1] <- sqrt(sigma2draw)
}
}
)[["elapsed"]]
## AR(1)-GARCH(1,1):
## allocate space:
realres[[3]] <- matrix(NA_real_, nrow = floor(draws/thinning), ncol = 3 + n + p)
colnames(realres[[3]]) <- c("a_0", "a_1", "b_1", paste("beta", 0:(p-1), sep="_"),
paste("sigma", 1:n, sep="_"))
## some auxiliary functions:
loglik <- function(y, s2) {
sum(dnorm(y, 0, sqrt(s2), log=TRUE))
}
s2s <- function(y, para, y20, s20) {
s2 <- rep(NA_real_, length(y))
s2[1] <- para[1] + para[2]*y20 + para[3]*s20
for (i in 2:length(y)) s2[i] <- para[1] + para[2]*y[i-1]^2 + para[3]*s2[i-1]
s2
}
## This is a very simple (and slow) random-walk MH sampler
## which can most certainly be improved...
## initial values:
# crude approximation of residual variance:
tmplm <- lm.fit(X, y)
s20 <- var(resid(tmplm))
y20 <- 0
## starting values:
## For the starting values for the parameters,
## use ML results from fGarch::garchFit on the
## residuals
# tmpfit <- fGarch::garchFit(~garch(1,1), resid(tmplm), trace=FALSE)
# startvals <- fGarch::coef(tmpfit)[c("omega", "alpha1", "beta1")]
startvals <- c(1.532868e-07, 3.246417e-02, 9.644149e-01)
## Random walk MH algorithm needs manual tuning for good mixing.
## Gelman et al. (1996) suggest an acceptance rate of 31.6%
## for spherical 3d multivariate normal (which we don't have :))
## These values seem to work OK:
# mhtuning <- tmpfit@fit$se.coef[c("omega", "alpha1", "beta1")] / 2
mhtuning <- c(3.463851e-08, 2.226946e-03, 2.354829e-03)
betadraw <- as.numeric(coef(tmplm))
paradraw <- startvals
paraprop <- rep(NA_real_, 3)
s2draw <- rep(.1, length(y))
s2prop <- rep(NA_real_, length(y))
accepts <- 0
## sampler:
tim3 <- system.time(
for (i in -(burnin-1):draws) {
if (i%%1000 == 0) cat("Iteration", i, "done.\n")
# draw volatilities and GARCH-parameters:
ytilde <- y - X %*% betadraw # regression residuals
paraprop <- rnorm(3, paradraw, mhtuning)
s2prop <- s2s(ytilde, paraprop, y20, s20)
if (all(s2prop > 0)) { # s2 needs to be positive, otherwise reject
logR <- loglik(ytilde, s2prop) - loglik(ytilde, s2draw) # log acceptance rate
if (log(runif(1)) < logR) {
paradraw <- paraprop
s2draw <- s2prop
if (i > 0) accepts <- accepts + 1
}
}
# draw the betas:
normalizer <- 1/sqrt(s2draw)
Xnew <- X * normalizer
ynew <- y * normalizer
Sigma <- solve(crossprod(Xnew) + B0inv)
mu <- Sigma %*% (crossprod(Xnew, ynew) + B0inv %*% b0)
betadraw <- as.numeric(mvtnorm::rmvnorm(1, mu, Sigma))
# store the results:
if (i > 0 & i %% thinning == 0) {
realres[[3]][i/thinning, paras] <- paradraw
realres[[3]][i/thinning, latents] <- sqrt(s2draw)
realres[[3]][i/thinning, betas] <- betadraw
}
}
)[["elapsed"]]
#Standardized residuals for model checking:
#homoskedastic
predmeans2 <- tcrossprod(X, realres[[2]][,1:2])
standresidmeans2 <- rowMeans((y - predmeans2)/realres[[2]][,3])
#heteroskedastic SV
predmeans <- tcrossprod(X, realres[[1]][,4:5])
standresidmeans <- rowMeans((y - predmeans)/exp(t(realres[[1]][,6:ncol(realres[[1]])])/2))
#heteroskedastic GARCH
predmeans3 <- tcrossprod(X, realres[[3]][,4:5])
standresidmeans3 <- rowMeans((y - predmeans3)/t(realres[[3]][,6:ncol(realres[[3]])]))
realresselection <- vector("list", 3)
realresselection[[1]] <- realres[[1]][,c('beta_0', 'beta_1')]
realresselection[[2]] <- realres[[2]][,c('beta_0', 'beta_1')]
realresselection[[3]] <- realres[[3]][,c('beta_0', 'beta_1')]
save(realresselection, standresidmeans, standresidmeans2, standresidmeans3, file = 'vignette_sampling_realworld.RData')
}
## ----betapost, fig.width=8.7, fig.height=7--------------------------
smootherfactor <- 1.8
par(mar = c(2.9, 2.9, 2.7, .5), mgp = c(1.8,.6,0), tcl = -.4)
layout(matrix(c(1,2,3,3), byrow=TRUE, nrow=2))
plot(density(realresselection[[1]][,"beta_0"], bw="SJ", adjust=smootherfactor),
main=bquote(paste("p(",beta[0],"|",bold(y),")", sep="")),
xlab="", ylab="", xlim=c(-.0005, .001))
lines(density(realresselection[[3]][,"beta_0"], bw="SJ", adjust=smootherfactor), lty=2, col=4)
lines(density(realresselection[[2]][,"beta_0"], bw="SJ", adjust=smootherfactor), lty=3, col=2)
#points(ols$coefficients[1],0, col=2)
#lines(confint(ols)[1,], rep(0,2), col=2)
legend("topleft", c("SV", "GARCH", "homosked."), col=c(1,4,2), lty=1:3)
plot(density(realresselection[[1]][,"beta_1"], bw="SJ", adjust=smootherfactor),
main=bquote(paste("p(",beta[1],"|",bold(y),")", sep="")),
xlab="", ylab="", xlim=c(0.9965, 1.0022))
lines(density(realresselection[[3]][,"beta_1"], bw="SJ", adjust=smootherfactor), lty=2, col=4)
lines(density(realresselection[[2]][,"beta_1"], bw="SJ", adjust=smootherfactor), lty=3, col=2)
#points(ols$coefficients[2],0, col=2)
#lines(confint(ols)[2,], rep(0,2), col=2)
legend("topright", c("SV", "GARCH", "homosked."), col=c(1,4,2), lty=1:3)
showevery <- 2L
plotorder <- sample.int(3*(nrow(realresselection[[1]]) / showevery))
cols <- rep(c("#000000bb", "#0000ff88", "#ff000088"),
each = nrow(realresselection[[1]]) / showevery)[plotorder]
pchs <- rep(1:3, each=nrow(realresselection[[1]]) / showevery)[plotorder]
myseq <- seq(1, nrow(realresselection[[1]]), by = showevery)
beta0 <- c(realresselection[[1]][myseq,"beta_0"],
realresselection[[3]][myseq,"beta_0"],
realresselection[[2]][myseq,"beta_0"])[plotorder]
beta1 <- c(realresselection[[1]][myseq,"beta_1"],
realresselection[[3]][myseq,"beta_1"],
realresselection[[2]][myseq,"beta_1"])[plotorder]
plot(beta0, beta1, col=cols, pch=pchs,
xlab=bquote(paste("p(",beta[0],"|", bold(y), ")")),
ylab=bquote(paste("p(",beta[1],"|", bold(y), ")")),
main="Scatterplot of posterior draws")
legend("topright", c("SV", "GARCH", "homosked."), col=c(1,4,2), pch=1:3)
## ----qqplot, fig.width=8, fig.height=7------------------------------
par(mfrow=c(3,2), mar = c(3.1, 3.3, 2.0, .5), mgp = c(1.7,.5,0),
tcl = -.4)
plot(standresidmeans2, main="Residual scatterplot (homoskedastic errors)",
xlab="Time", ylab="Standardized residuals")
qqplot(qnorm(ppoints(length(standresidmeans2)), mean=0, sd=1),
standresidmeans2, main="Residual Q-Q plot (homoskedastic errors)",
xlab = "Theoretical N(0,1)-quantiles", ylab = "Empirical Quantiles")
abline(0,1)
plot(standresidmeans3, main="Residual scatterplot (GARCH errors)",
xlab="Time", ylab="Standardized residuals")
qqplot(qnorm(ppoints(length(standresidmeans3)), mean=0, sd=1),
standresidmeans3, main="Residual Q-Q plot (GARCH errors)",
xlab = "Theoretical N(0,1)-quantiles", ylab = "Empirical Quantiles")
abline(0,1)
plot(standresidmeans, main="Residual scatterplot (SV errors)",
xlab="Time", ylab="Standardized residuals")
qqplot(qnorm(ppoints(length(standresidmeans)), mean=0, sd=1),
standresidmeans, main="Residual Q-Q plot (SV errors)",
xlab = "Theoretical N(0,1)-quantiles", ylab = "Empirical Quantiles")
abline(0,1)
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.