# Empirical Bayes

library(bvhar)

# Hyperparameters in Normal-inverse-Wishart Priors

In this vignette, we discuss hyperparameters of Normal-inverse-Wishart priors for BVAR, BVHAR-S, and BVHAR-L.

## BVAR

From Litterman (1986) and Bańbura et al. (2010), there are $$\sigma_j$$ (sigma), $$\lambda$$ (lambda), and $$\delta_j$$ (delta) in BVAR.

• $$\sigma_i^2 / \sigma_j^2$$ in Minnesota moments explain the data scales.
• $$\delta_j$$’s are related to the belief to random walk.
• If $$\forall j \; \delta_j = 1$$, it is random walk prior (Litterman (1986) setting).
• If $$\forall j \; \delta_j = 0$$, it is white noise.
• $$\lambda$$ controls the overall tightness of the prior around these two prior beliefs.

In addition, there is small number $$\epsilon$$ (eps) for matrix invertibility. This is set to be 1e-04 (default). You can change the number. Based on these properties, users should subjectively choose hyperparameter.

Consider four indices of CBOE ETF volatility index (etf_vix): Gold, crude oil, emerging markets, gold miners, and silver. We split train-test set (Test = last 20 points).

etf_split <-
etf_vix %>%
dplyr::select(GVZCLS, OVXCLS, VXEEMCLS, VXGDXCLS, VXSLVCLS) %>%
divide_ts(20)
# split---------------------
etf_train <- etf_split$train etf_test <- etf_split$test
dim_data <- ncol(etf_train)

The following is BVAR Minnesota prior specification.

• sigma: standard error
• lambda: not that small number because small dimension
• delta: Litterman’s setting
sig <- apply(etf_train, 2, sd)
lam <- .2
del <- rep(1, dim_data)
# bvharspec------------------
(bvar_spec <- set_bvar(
sigma = sig,
lambda = lam,
delta = del,
eps = 1e-04
))
#> Model Specification for BVAR
#>
#> Parameters: Coefficent matrice and Covariance matrix
#> Prior: Minnesota
#> # Type '?bvar_minnesota' in the console for some help.
#> ========================================================
#>
#> Setting for 'sigma':
#>   GVZCLS    OVXCLS  VXEEMCLS  VXGDXCLS  VXSLVCLS
#>     3.77     10.63      4.39      7.45      5.99
#>
#> Setting for 'lambda':
#> [1]  0.2
#>
#> Setting for 'delta':
#> [1]  1  1  1  1  1
#>
#> Setting for 'eps':
#> [1]  1e-04

There are one more parameter, order $$p$$. Set this $$p = 3$$.

bvar_cand <- bvar_minnesota(
y = etf_train,
p = 3,
bayes_spec = bvar_spec,
include_mean = TRUE
)

## BVHAR-S

BVHAR-S (Kim et al. (n.d.)) has the same set of hyperparameters with BVAR.

(bvhar_short_spec <- set_bvhar(
sigma = sig,
lambda = lam,
delta = del,
eps = 1e-04
))
#> Model Specification for BVHAR
#>
#> Parameters: Coefficent matrice and Covariance matrix
#> Prior: MN_VAR
#> # Type '?bvhar_minnesota' in the console for some help.
#> ========================================================
#>
#> Setting for 'sigma':
#>   GVZCLS    OVXCLS  VXEEMCLS  VXGDXCLS  VXSLVCLS
#>     3.77     10.63      4.39      7.45      5.99
#>
#> Setting for 'lambda':
#> [1]  0.2
#>
#> Setting for 'delta':
#> [1]  1  1  1  1  1
#>
#> Setting for 'eps':
#> [1]  1e-04
bvhar_short_cand <- bvhar_minnesota(
y = etf_train,
har = c(5, 22),
bayes_spec = bvhar_short_spec,
include_mean = TRUE
)

## BVHAR-L

While above two priors set prior mean zero for farther coefficient (i.e. $$A_2, \ldots, A_p$$ or $$\Phi^{(w)}, \Phi^{(m)}$$), now these weekly and monthly coefficients are also have nonzero prior mean. So instead of $$\delta_j$$, there are

• $$d_j$$ for daily term.
• $$w_j$$ for weekly term.
• $$m_j$$ for monthly term.
dayj <- rep(.8, dim_data)
weekj <- rep(.2, dim_data)
monthj <- rep(.1, dim_data)
# bvharspec------------------
bvhar_long_spec <- set_weight_bvhar(
sigma = sig,
lambda = lam,
eps = 1e-04,
daily = dayj,
weekly = weekj,
monthly = monthj
)
bvhar_long_cand <- bvhar_minnesota(
y = etf_train,
har = c(5, 22),
bayes_spec = bvhar_long_spec,
include_mean = TRUE
)

# Hyperparameter Selection

• Giannone et al. (2015) provides the closed form of marginal likelihood for BVAR Minnesota prior.
• Based on this calculation, Kim et al. (n.d.) gives the closed form for BVHAR-S and BVHAR-L.

This package provides hyperparameter selection function with empirical bayes (Giannone et al. (2015)). Implementation is simple. Provide above bvharspec as initial values to stats::optim() structure. You can either use

• Individual choose_bvar() or choose_bvhar() function
• or one integrated choose_bayes() function

## Individual Functions

You can parallelize the computation with optimParallel::optimParallel() by specifying parallel = list() (if you leave it as empty list, the function execute stats::optim()). Register cluster, and pass cl = cl. For the details of other two, please see the documentation of the optimParallel::optimParallel().

• cl: Register cluster
• forward: Recommendation - FALSE
• loginfo: Recommendation - FALSE
cl <- parallel::makeCluster(2)

### BVAR

We first try BVAR. choose_bvar(bayes_spec, lower = .01, upper = 10, eps = 1e-04, y, p, include_mean = TRUE, parallel = list()) chooses BVAR hyperparameter. Observe that it fixes the order p (of course eps, either).

By default, it apply .01 to lower bound and 10 to upper bound. When using the function, setting lower and upper bounds can be quite tricky. It’s because the bound should be expressed by vector as in the stats::optim() function. See the following code how to specify the bound.

(bvar_optim <- choose_bvar(
bayes_spec = bvar_spec,
lower = c(
rep(1, dim_data), # sigma
1e-2, # lambda
rep(1e-2, dim_data) # delta
),
upper = c(
rep(15, dim_data), # sigma
Inf, # lambda
rep(1, dim_data) # delta
),
eps = 1e-04,
y = etf_train,
p = 3,
include_mean = TRUE,
parallel = list(cl = cl, forward = FALSE, loginfo = FALSE)
))
#> Model Specification for BVAR
#>
#> Parameters: Coefficent matrice and Covariance matrix
#> Prior: Minnesota
#> # Type '?bvar_minnesota' in the console for some help.
#> ========================================================
#>
#> Setting for 'sigma':
#>   GVZCLS    OVXCLS  VXEEMCLS  VXGDXCLS  VXSLVCLS
#>     2.50      3.91      3.54      4.51      3.78
#>
#> Setting for 'lambda':
#>
#> 0.579
#>
#> Setting for 'delta':
#>
#> 1  1  1  1  1
#>
#> Setting for 'eps':
#> [1]  1e-04

The result is a class named bvharemp.

class(bvar_optim)
#> [1] "bvharemp"
names(bvar_optim)
#> [1] "par"         "value"       "counts"      "convergence" "message"
#> [6] "spec"        "fit"         "ml"

It has the chosen specification (bvharspec), fit (bvarmn), and its marginal likelihood value (ml).

### BVHAR-L

Using the same choose_bvhar() function, you should change bayes_spec and the length of bounds vectors:

(bvhar_long_optim <- choose_bvhar(
bayes_spec = bvhar_long_spec,
lower = c(
rep(1, dim_data), # sigma
1e-2, # lambda
rep(1e-2, dim_data), # daily
rep(1e-2, dim_data), # weekly
rep(1e-2, dim_data) # monthly
),
upper = c(
rep(15, dim_data), # sigma
Inf, # lambda
rep(1, dim_data), # daily
rep(1, dim_data), # weekly
rep(1, dim_data) # monthly
),
eps = 1e-04,
y = etf_train,
har = c(5, 22),
include_mean = TRUE,
parallel = list(cl = cl, forward = FALSE, loginfo = FALSE)
))
#> Model Specification for BVHAR
#>
#> Parameters: Coefficent matrice and Covariance matrix
#> Prior: MN_VHAR
#> # Type '?bvhar_minnesota' in the console for some help.
#> ========================================================
#>
#> Setting for 'sigma':
#>   GVZCLS    OVXCLS  VXEEMCLS  VXGDXCLS  VXSLVCLS
#>     2.53      3.94      3.55      5.06      4.26
#>
#> Setting for 'lambda':
#>
#> 0.519
#>
#> Setting for 'eps':
#> [1]  1e-04
#>
#> Setting for 'daily':
#>
#> 1  1  1  1  1
#>
#> Setting for 'weekly':
#>
#> 0.316  0.226  0.234  0.806  0.542
#>
#> Setting for 'monthly':
#>
#> 0.2184  0.1444  0.0995  0.2681  0.2399
class(bvhar_long_optim)
#> [1] "bvharemp"
names(bvhar_long_optim)
#> [1] "par"         "value"       "counts"      "convergence" "message"
#> [6] "spec"        "fit"         "ml"

Final fit:

bvhar_long_final <- bvhar_long_optim\$fit

## Integrated Function

As mentioned, choose_bvar() and choose_bvhar() are still providing difficult ways of bounding methods. So, there are another function: choose_bayes(). You can set Empirical Bayes bound using bound_bvhar(init_spec, lower_spec, upper_spec).

# lower bound----------------
bvar_lower <- set_bvar(
sigma = rep(1, dim_data),
lambda = 1e-2,
delta = rep(1e-2, dim_data)
)
# upper bound---------------
bvar_upper <- set_bvar(
sigma = rep(15, dim_data),
lambda = Inf,
delta = rep(1, dim_data)
)
# bound--------------------
(bvar_bound <- bound_bvhar(
init_spec = bvar_spec,
lower_spec = bvar_lower,
upper_spec = bvar_upper
))
#> Lower bound specification for BVAR optimization (L-BFGS-B):
#> ========================================================
#> ========================================================
#> Model Specification for BVAR
#>
#> Parameters: Coefficent matrice and Covariance matrix
#> Prior: Minnesota
#> # Type '?bvar_minnesota' in the console for some help.
#> ========================================================
#>
#> Setting for 'sigma':
#> [1]  1  1  1  1  1
#>
#> Setting for 'lambda':
#> [1]  0.01
#>
#> Setting for 'delta':
#> [1]  0.01  0.01  0.01  0.01  0.01
#>
#> Setting for 'eps':
#> [1]  1e-04
#>
#>
#>
#> Upper bound specification for BVAR optimization (L-BFGS-B):
#> ========================================================
#> ========================================================
#> Model Specification for BVAR
#>
#> Parameters: Coefficent matrice and Covariance matrix
#> Prior: Minnesota
#> # Type '?bvar_minnesota' in the console for some help.
#> ========================================================
#>
#> Setting for 'sigma':
#> [1]  15  15  15  15  15
#>
#> Setting for 'lambda':
#> [1]  Inf
#>
#> Setting for 'delta':
#> [1]  1  1  1  1  1
#>
#> Setting for 'eps':
#> [1]  1e-04
class(bvar_bound)
#> [1] "boundbvharemp"
names(bvar_bound)
#> [1] "spec"  "lower" "upper"

Based on this boundbvharemp, we can use choose_bayes() function. This function implements choose_bvar() or choose_bvhar() given inputs.

• bayes_bound: boundbvharemp object.
• order: p of BVAR or har of BVHAR.
• The other options are the same.
(bvar_optim_v2 <- choose_bayes(
bayes_bound = bvar_bound,
eps = 1e-04,
y = etf_train,
order = 3,
include_mean = TRUE,
parallel = list(cl = cl, forward = FALSE, loginfo = FALSE)
))
#> Model Specification for BVAR
#>
#> Parameters: Coefficent matrice and Covariance matrix
#> Prior: Minnesota
#> # Type '?bvar_minnesota' in the console for some help.
#> ========================================================
#>
#> Setting for 'sigma':
#>   GVZCLS    OVXCLS  VXEEMCLS  VXGDXCLS  VXSLVCLS
#>     2.50      3.91      3.54      4.51      3.78
#>
#> Setting for 'lambda':
#>
#> 0.579
#>
#> Setting for 'delta':
#>
#> 1  1  1  1  1
#>
#> Setting for 'eps':
#> [1]  1e-04

Do not forget shut down the cluster cl.

parallel::stopCluster(cl)