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()
orchoose_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
: Registercluster
-
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.90 3.54 4.50 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.05
#>
#> Setting for 'lambda':
#>
#> 0.551
#>
#> 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.05 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.225 0.234 0.805 0.542
#>
#> Setting for 'monthly':
#>
#> 0.2180 0.1438 0.0998 0.2674 0.2408
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
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 orhar
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.90 3.54 4.50 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)