ier {bridgesampling} | R Documentation |
Standardized International Exchange Rate Changes from 1975 to 1986
Description
This data set contains the changes in monthly international exchange rates for pounds sterling from January 1975 to December 1986 obtained from West and Harrison (1997, pp. 612-615). Currencies tracked are US Dollar (column us_dollar
), Canadian Dollar (column canadian_dollar
), Japanese Yen (column yen
), French Franc (column franc
), Italian Lira (column lira
), and the (West) German Mark (column mark
). Each series has been standardized with respect to its sample mean and standard deviation.
Usage
ier
Format
A matrix with 143 rows and 6 columns.
Source
West, M., Harrison, J. (1997). Bayesian forecasting and dynamic models (2nd ed.). Springer-Verlag, New York.
Lopes, H. F., West, M. (2004). Bayesian model assessment in factor analysis. Statistica Sinica, 14, 41-67. https://www.jstor.org/stable/24307179
Examples
## Not run:
################################################################################
# BAYESIAN FACTOR ANALYSIS (AS PROPOSED BY LOPES & WEST, 2004)
################################################################################
library(bridgesampling)
library(rstan)
cores <- 4
options(mc.cores = cores)
data("ier")
#-------------------------------------------------------------------------------
# plot data
#-------------------------------------------------------------------------------
currency <- colnames(ier)
label <- c("US Dollar", "Canadian Dollar", "Yen", "Franc", "Lira", "Mark")
op <- par(mfrow = c(3, 2), mar = c(6, 6, 3, 3))
for (i in seq_along(currency)) {
plot(ier[,currency[i]], type = "l", col = "darkblue", axes = FALSE,
ylim = c(-4, 4), ylab = "", xlab = "", lwd = 2)
axis(1, at = 0:12*12, labels = 1975:1987, cex.axis = 1.7)
axis(2, at = pretty(c(-4, 4)), las = 1, cex.axis = 1.7)
mtext("Year", 1, cex = 1.5, line = 3.2)
mtext("Exchange Rate Changes", 2, cex = 1.4, line = 3.2)
mtext(label[i], 3, cex = 1.6, line = .1)
}
par(op)
#-------------------------------------------------------------------------------
# stan model
#-------------------------------------------------------------------------------
model_code <-
"data {
int<lower=1> T; // number of observations
int<lower=1> m; // number of variables
int<lower=1> k; // number of factors
matrix[T,m] Y; // data matrix
}
transformed data {
int<lower = 1> r;
vector[m] zeros;
r = m * k - k * (k - 1) / 2; // number of non-zero factor loadings
zeros = rep_vector(0.0, m);
}
parameters {
real beta_lower[r - k]; // lower-diagonal elements of beta
real<lower = 0> beta_diag [k]; // diagonal elements of beta
vector<lower = 0>[m] sigma2; // residual variances
}
transformed parameters {
matrix[m,k] beta;
cov_matrix[m] Omega;
// construct lower-triangular factor loadings matrix
{
int index_lower = 1;
for (j in 1:k) {
for (i in 1:m) {
if (i == j) {
beta[j,j] = beta_diag[j];
} else if (i >= j) {
beta[i,j] = beta_lower[index_lower];
index_lower = index_lower + 1;
} else {
beta[i,j] = 0.0;
}
}
}
}
Omega = beta * beta' + diag_matrix(sigma2);
}
model {
// priors
target += normal_lpdf(beta_diag | 0, 1) - k * normal_lccdf(0 | 0, 1);
target += normal_lpdf(beta_lower | 0, 1);
target += inv_gamma_lpdf(sigma2 | 2.2 / 2.0, 0.1 / 2.0);
// likelihood
for(t in 1:T) {
target += multi_normal_lpdf(Y[t] | zeros, Omega);
}
}"
# compile model
model <- stan_model(model_code = model_code)
#-------------------------------------------------------------------------------
# fit models and compute log marginal likelihoods
#-------------------------------------------------------------------------------
# function for generating starting values
init_fun <- function(nchains, k, m) {
r <- m * k - k * (k - 1) / 2
out <- vector("list", nchains)
for (i in seq_len(nchains)) {
beta_lower <- array(runif(r - k, 0.05, 1), dim = r - k)
beta_diag <- array(runif(k, .05, 1), dim = k)
sigma2 <- array(runif(m, .05, 1.5), dim = m)
out[[i]] <- list(beta_lower = beta_lower,
beta_diag = beta_diag,
sigma2 = sigma2)
}
return(out)
}
set.seed(1)
stanfit <- bridge <- vector("list", 3)
for (k in 1:3) {
stanfit[[k]] <- sampling(model,
data = list(Y = ier, T = nrow(ier),
m = ncol(ier), k = k),
iter = 11000, warmup = 1000, chains = 4,
init = init_fun(nchains = 4, k = k, m = ncol(ier)),
cores = cores, seed = 1)
bridge[[k]] <- bridge_sampler(stanfit[[k]], method = "warp3",
repetitions = 10, cores = cores)
}
# example output
summary(bridge[[2]])
#-------------------------------------------------------------------------------
# compute posterior model probabilities
#-------------------------------------------------------------------------------
pp <- post_prob(bridge[[1]], bridge[[2]], bridge[[3]],
model_names = c("k = 1", "k = 2", "k = 3"))
pp
op <- par(mar = c(6, 6, 3, 3))
boxplot(pp, axes = FALSE,
ylim = c(0, 1), ylab = "",
xlab = "")
axis(1, at = 1:3, labels = colnames(pp), cex.axis = 1.7)
axis(2, cex.axis = 1.1)
mtext("Posterior Model Probability", 2, cex = 1.5, line = 3.2)
mtext("Number of Factors", 1, cex = 1.4, line = 3.2)
par(op)
## End(Not run)