# Gibbs sampler for vector autoregressive model. Obtain samples of the posterior of a Bayesian VAR model of fixed order. An independent Normal-Inverse-Wishart prior is employed. ```r gibbs_var( data, ar.order, Ntotal, burnin, thin = 1, print_interval = 500, full_lik = F, beta.mu = rep(0, ar.order * ncol(data)^2), beta.Sigma = 10000 * diag(ar.order * ncol(data)^2), Sigma.S = 1e-04 * diag(ncol(data)), Sigma.nu = 1e-04 ) ``` ## Arguments - `data`: numeric matrix; NA values are interpreted as missing values and treated as random - `ar.order`: order of the autoregressive model (integer >= 0) - `Ntotal`: total number of iterations to run the Markov chain - `burnin`: number of initial iterations to be discarded - `thin`: thinning number (postprocessing) - `print_interval`: Number of iterations, after which a status is printed to console - `full_lik`: logical; if TRUE, the full likelihood for all observations is used; if FALSE, the partial likelihood for the last n-p observations - `beta.mu`: prior mean of beta vector (normal) - `beta.Sigma`: prior covariance matrix of beta vector - `Sigma.S`: prior parameter for the innovation covariance matrix, symmetric positive definite matrix - `Sigma.nu`: prior parameter for the innovation covariance matrix, nonnegative real number ## Returns list containing the following fields: - **beta**: matrix containing traces of the VAR parameter vector beta - **Sigma**: trace of innovation covariance Sigma - **psd.median,psd.mean**: psd estimates: (pointwise, componentwise) posterior median and mean - **psd.p05,psd.p95**: pointwise credibility interval - **psd.u05,psd.u95**: uniform credibility interval, see (6.5) in Meier (2018) - **lpost**: trace of log posterior ## Details See Section 2.2.3 in Koop and Korobilis (2010) or Section 6.2 in Meier (2018) for further details ## Examples ```r ## Not run: ## ## Example 1: Fit a VAR(p) model to SOI/Recruitment series: ## # Use this variable to set the VAR model order p <- 5 data <- cbind(as.numeric(astsa::soi-mean(astsa::soi)), as.numeric(astsa::rec-mean(astsa::rec)) / 50) data <- apply(data, 2, function(x) x-mean(x)) # If you run the example be aware that this may take several minutes print("example may take some time to run") mcmc <- gibbs_var(data=data, ar.order=p, Ntotal=10000, burnin=4000, thin=2) # Plot spectral estimate, credible regions and periodogram on log-scale plot(mcmc, log=T) ## ## Example 2: Fit a VAR(p) model to VMA(1) data ## # Use this variable to set the VAR model order p <- 5 n <- 256 ma <- rbind(c(-0.75, 0.5), c(0.5, 0.75)) Sigma <- rbind(c(1, 0.5), c(0.5, 1)) data <- sim_varma(model=list(ma=ma), n=n, d=2) data <- apply(data, 2, function(x) x-mean(x)) # If you run the example be aware that this may take several minutes print("example may take some time to run") mcmc <- gibbs_var(data=data, ar.order=p, Ntotal=10000, burnin=4000, thin=2) # Plot spectral estimate, credible regions and periodogram on log-scale plot(mcmc, log=T) ## End(Not run) ``` ## References G. Koop and D. Korobilis (2010) **Bayesian Multivariate Time Series Methods for Empirical Macroeconomics** Foundations and Trends in Econometrics A. Meier (2018) **A Matrix Gamma Process and Applications to Bayesian Analysis of Multivariate Time Series** PhD thesis, OvGU Magdeburg

gibbs_var function