## ar.R -- to accompany Question 12 of the second example sheet ## for Part III/Mphil Monte Carlo Inference, Lent 2008, ## Mathematical Tripos ## ## by Robert B. Gramacy, Statistical Laboratory ## University of Cambridge ## get a function for drawing from multivariate normal distributions if(! suppressWarnings(require(mvtnorm, quietly=TRUE))) { source("rmultnorm.R") rmvnorm <- rmultnorm } ## a.gibbs: ## ## single gibbs sample for the AR coefficients a a.gibbs <- function(x, Xk, se2, mu.a, sigma2.a) { ## find k and n k <- ncol(Xk) if(is.null(k)) k <- 1 ## find C, the k-MVN covariance matrix Ci <- (t(Xk) %*% Xk)/se2 + diag(1/sigma2.a,nrow=k) C <- solve(Ci) ## find mu, the k-MVN mean mu <- C %*% ((t(Xk) %*% x)/se2 - mu.a/sigma2.a) # draw ak from the k-MVN ak <- rmvnorm(1, mu, C) return(ak) } ## se2.gibbs: ## ## single gibbs sample for the AR variance parameter se2 se2.gibbs <- function(x, a, Xk, alpha, beta) { ## sum of squared differences of length k=length(epsilon) epsilon <- x - Xk %*% a nmkmax <- length(epsilon) ete <- sum(epsilon^2) ## inverse gamma parameters a <- (alpha + nmkmax)/2 b <- (ete + beta)/2 ## draw from an inverse gamma distribution return(1/rgamma(1, shape=a, scale=1/b)) } ## ar.gibbs: ## ## Gibbs Sampling for Autogregressive models of order ## k = ncol(Xk) -- i.e., specified in matrix form ar.gibbs <- function(T, ## total number of MCMC rounds x, Xk, ## input data and model order mu.a=NULL, sigma2.a=10^2, ## mu priors alpha=2*(10^(-3)), beta=2*(10^(-3)), ## sigma2 priors quiet=TRUE) { ## infer k k <- ncol(Xk) if(is.null(k)) k <- 1 ## make mu.a into a vector if(is.null(mu.a)) mu.a <- rep(0, k) else if(length(mu.a) != k) stop("bad length(mu.a)") ## allocate collections of parameters to return se2 <- rep(1, T) a <- matrix(0, nrow=T, ncol=k) ## each round of Gibbs Sampling for(t in 2:T) { ## gibbs draw for se2 se2[t] <- se2.gibbs(x, a[t-1,], Xk, alpha, beta) ## gibbs draw for a a[t,] <- as.vector(a.gibbs(x, Xk, se2[t], mu.a, sigma2.a)) ## print progress meter if(!quiet && t%%1000 == 0) cat(paste("t=", t, "\n", sep="")) } ##return return(data.frame(a=a, se2=se2)) } ## make.XKs ## ## construct the Xk matrices for each k=1,...,kmax from the ## time series vector x. make.XKs <- function(x, kmax) { n <- length(x) XKs <- list() XKs[[1]] <- matrix(x[kmax:(n-1)], ncol=1) for(k in 2:kmax) { XKs[[k]] <- as.matrix(cbind(XKs[[k-1]], x[(kmax-k+1):(n-k)])) } return(XKs) }