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).

bvar_final <- bvar_optim$fit

BVHAR-S

choose_bvhar(bayes_spec, lower = .01, upper = 10, eps = 1e-04, y, har = c(5, 22), include_mean = TRUE, parallel = list()) choses hyperparameter set of BVHAR-S or BVHAR-L given bayes_spec. If it is set_bvhar(), it performs BVHAR-S. It it is set_weight_bvhar(), it performs BVHAR-L.

(bvhar_short_optim <- choose_bvhar(
  bayes_spec = bvhar_short_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,
  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_VAR
#> # Type '?bvhar_minnesota' in the console for some help.
#> ========================================================
#> 
#> Setting for 'sigma':
#>   GVZCLS    OVXCLS  VXEEMCLS  VXGDXCLS  VXSLVCLS  
#>     2.54      4.07      3.67      4.43      4.06  
#> 
#> Setting for 'lambda':
#>        
#> 0.552  
#> 
#> Setting for 'delta':
#>                
#> 1  1  1  1  1  
#> 
#> Setting for 'eps':
#> [1]  1e-04

The structure of the result is the same.

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

bvhar_short_final <- bvhar_short_optim$fit

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)