Normal-inverse-Wishart Matrix

We provide functions to generate matrix-variate Normal and inverse-Wishart.

sim_mnormal(num_sim, mu, sig): num_sim of \(\mathbf{X}_i \stackrel{iid}{\sim} N(\boldsymbol{\mu}, \Sigma)\).
sim_matgaussian(mat_mean, mat_scale_u, mat_scale_v): One \(X_{m \times n} \sim MN(M_{m \times n}, U_{m \times m}, V_{n \times n})\) which means that \(vec(X) \sim N(vec(M), V \otimes U)\).
sim_iw(mat_scale, shape): One \(\Sigma \sim IW(\Psi, \nu)\).
sim_mniw(num_sim, mat_mean, mat_scale_u, mat_scale, shape): num_sim of \((X_i, \Sigma_i) \stackrel{iid}{\sim} MNIW(M, U, V, \nu)\).

Multivariate Normal generation gives num_sim x dim matrix. For example, generating 3 vector from Normal(\(\boldsymbol{\mu} = \mathbf{0}_2\), \(\Sigma = diag(\mathbf{1}_2)\)):

sim_mnormal(3, rep(0, 2), diag(2))
#>        [,1]   [,2]
#> [1,] -0.626  0.184
#> [2,] -0.836  1.595
#> [3,]  0.330 -0.820

The output of sim_matgaussian() is a matrix.

sim_matgaussian(matrix(1:20, nrow = 4), diag(4), diag(5), FALSE)
#>      [,1] [,2]  [,3] [,4] [,5]
#> [1,] 1.49 5.74  9.58 12.7 18.5
#> [2,] 2.39 5.38  7.79 15.1 18.0
#> [3,] 2.98 7.94 11.82 15.6 19.9
#> [4,] 4.78 8.07 10.01 16.6 19.9

When generating IW, violating \(\nu > dim - 1\) gives error. But we ignore \(\nu > dim + 1\) (condition for mean existence) in this function. Nonetheless, we recommend you to keep \(\nu > dim + 1\) condition. As mentioned, it guarantees the existence of the mean.

sim_iw(diag(5), 7)
#>         [,1]    [,2]    [,3]    [,4]    [,5]
#> [1,]  0.1827  0.0894 -0.0411 -0.0924 -0.1305
#> [2,]  0.0894  0.6110  0.0860 -0.3754 -0.1015
#> [3,] -0.0411  0.0860  0.2189 -0.1649 -0.0988
#> [4,] -0.0924 -0.3754 -0.1649  0.4577  0.2136
#> [5,] -0.1305 -0.1015 -0.0988  0.2136  0.7444

In case of sim_mniw(), it returns list with mn (stacked MN matrices) and iw (stacked IW matrices). Each mn and iw has draw lists.

sim_mniw(2, matrix(1:20, nrow = 4), diag(4), diag(5), 7, FALSE)
#> $mn
#> $mn[[1]]
#>      [,1] [,2]  [,3] [,4] [,5]
#> [1,] 1.19 4.96  8.03 12.7 14.6
#> [2,] 2.13 5.17 10.41 14.1 19.3
#> [3,] 3.01 7.15 10.29 14.9 18.1
#> [4,] 4.59 7.99 12.12 15.0 18.5
#> 
#> $mn[[2]]
#>       [,1] [,2]  [,3] [,4] [,5]
#> [1,] 0.323 5.07  8.24 13.1 17.4
#> [2,] 2.022 5.56  9.96 13.9 19.4
#> [3,] 3.380 7.24 11.34 15.5 17.9
#> [4,] 3.755 8.79 11.84 15.6 19.6
#> 
#> 
#> $iw
#> $iw[[1]]
#>         [,1]     [,2]     [,3]    [,4]   [,5]
#> [1,]  0.2887 -0.04393  0.11892 -0.2959  0.110
#> [2,] -0.0439  0.33357  0.00197  0.0106 -0.167
#> [3,]  0.1189  0.00197  1.30074 -0.0291  2.210
#> [4,] -0.2959  0.01061 -0.02913  0.5072  0.250
#> [5,]  0.1104 -0.16736  2.20957  0.2504  4.623
#> 
#> $iw[[2]]
#>          [,1]    [,2]     [,3]    [,4]   [,5]
#> [1,]  0.28118 -0.1397  0.00529  0.0283  0.038
#> [2,] -0.13965  0.2876  0.06156 -0.1575 -0.186
#> [3,]  0.00529  0.0616  0.31282 -0.1372 -0.299
#> [4,]  0.02831 -0.1575 -0.13724  0.3441  0.122
#> [5,]  0.03803 -0.1856 -0.29860  0.1217  0.738

This function has been defined for the next simulation functions.

Minnesota Prior

BVAR

Consider BVAR Minnesota prior setting,

\[A \sim MN(A_0, \Omega_0, \Sigma_e)\]

\[\Sigma_e \sim IW(S_0, \alpha_0)\]

From Litterman (1986) and Bańbura et al. (2010)
Each \(A_0, \Omega_0, S_0, \alpha_0\) is defined by adding dummy observations
- build_xdummy()
- build_ydummy()
sigma: Vector \(\sigma_1, \ldots, \sigma_m\)
- \(\Sigma_e = diag(\sigma_1^2, \ldots, \sigma_m^2)\)
- \(\sigma_i^2 / \sigma_j^2\): different scale and variability of the data
lambda
- Controls the overall tightness of the prior distribution around the RW or WN
- Governs the relative importance of the prior beliefs w.r.t. the information contained in the data
  - If \(\lambda = 0\), then posterior = prior and the data do not influence the estimates.
  - If \(\lambda = \infty\), then posterior expectations = OLS.
- Choose in relation to the size of the system (Bańbura et al. (2010))
  - As m increases, \(\lambda\) should be smaller to avoid overfitting (De Mol et al. (2008))
delta: Persistence
- Litterman (1986) originally sets high persistence \(\delta_i = 1\)
- For Non-stationary variables: random walk prior \(\delta_i = 1\)
- For stationary variables: white noise prior \(\delta_i = 0\)
eps: Very small number to make matrix invertible

bvar_lag <- 5
(spec_to_sim <- set_bvar(
  sigma = c(3.25, 11.1, 2.2, 6.8), # sigma vector
  lambda = .2, # lambda
  delta = rep(1, 4), # 4-dim delta vector
  eps = 1e-04 # very small number
))
#> Model Specification for BVAR
#> 
#> Parameters: Coefficent matrice and Covariance matrix
#> Prior: Minnesota
#> ========================================================
#> 
#> Setting for 'sigma':
#> [1]   3.25  11.10   2.20   6.80
#> 
#> Setting for 'lambda':
#> [1]  0.2
#> 
#> Setting for 'delta':
#> [1]  1  1  1  1
#> 
#> Setting for 'eps':
#> [1]  1e-04
#> 
#> Setting for 'hierarchical':
#> [1]  FALSE

sim_mncoef(p, bayes_spec, full = TRUE) can generate both \(A\) and \(\Sigma\) matrices.
In bayes_spec, only set_bvar() works.
If full = FALSE, \(\Sigma\) is not random. It is same as diag(sigma) from the bayes_spec.
full = TRUE is the default.

(sim_mncoef(bvar_lag, spec_to_sim))
#> $coefficients
#>            [,1]     [,2]      [,3]     [,4]
#>  [1,]  0.996255 -0.19084 -0.026363  0.07600
#>  [2,] -0.017158  1.07873 -0.033404 -0.00816
#>  [3,] -0.225869  0.31165  0.927705 -0.48148
#>  [4,] -0.002706 -0.26068  0.038566  0.84105
#>  [5,] -0.023064 -0.11717 -0.030881  0.35496
#>  [6,]  0.000253 -0.05523 -0.018351  0.03580
#>  [7,] -0.001364 -0.06563 -0.051615 -0.33330
#>  [8,] -0.025569  0.11938 -0.011720 -0.15248
#>  [9,]  0.062006 -0.07985 -0.001830  0.15440
#> [10,]  0.008609 -0.03437  0.016242  0.00866
#> [11,] -0.069804  0.11559  0.003275  0.30043
#> [12,]  0.001677  0.01789 -0.000582 -0.02224
#> [13,] -0.000864  0.05893  0.038960  0.06788
#> [14,]  0.008799 -0.04337  0.005557  0.02273
#> [15,] -0.053924  0.16120 -0.006870  0.10085
#> [16,] -0.009109  0.00284 -0.008880 -0.01168
#> [17,]  0.008001 -0.05816  0.012293 -0.05256
#> [18,] -0.006107  0.03645 -0.004491 -0.00687
#> [19,] -0.011204  0.05703  0.036651  0.13185
#> [20,] -0.003148 -0.05735  0.016794  0.04086
#> 
#> $covmat
#>        [,1]     [,2]    [,3]   [,4]
#> [1,]  2.615  -5.1044  0.1522   2.44
#> [2,] -5.104  32.2768 -0.0549 -12.04
#> [3,]  0.152  -0.0549  1.4573   0.31
#> [4,]  2.440 -12.0446  0.3101  30.83

BVHAR

sim_mnvhar_coef(bayes_spec, full = TRUE) generates BVHAR model setting:

\[\Phi \mid \Sigma_e \sim MN(M_0, \Omega_0, \Sigma_e)\]

\[\Sigma_e \sim IW(\Psi_0, \nu_0)\]

similar to BVAR, bayes_spec option wants bvharspec. But
- set_bvhar()
- set_weight_bvhar()
There is full = TRUE, too.

BVHAR-S

(bvhar_var_spec <- set_bvhar(
  sigma = c(1.2, 2.3), # sigma vector
  lambda = .2, # lambda
  delta = c(.3, 1), # 2-dim delta vector
  eps = 1e-04 # very small number
))
#> Model Specification for BVHAR
#> 
#> Parameters: Coefficent matrice and Covariance matrix
#> Prior: MN_VAR
#> ========================================================
#> 
#> Setting for 'sigma':
#> [1]  1.2  2.3
#> 
#> Setting for 'lambda':
#> [1]  0.2
#> 
#> Setting for 'delta':
#> [1]  0.3  1.0
#> 
#> Setting for 'eps':
#> [1]  1e-04
#> 
#> Setting for 'hierarchical':
#> [1]  FALSE

(sim_mnvhar_coef(bvhar_var_spec))
#> $coefficients
#>          [,1]    [,2]
#> [1,]  0.30349 -0.2954
#> [2,] -0.01177  1.2277
#> [3,]  0.05543 -0.2689
#> [4,]  0.00235  0.0574
#> [5,] -0.05121  0.1928
#> [6,] -0.00891  0.0245
#> 
#> $covmat
#>        [,1]   [,2]
#> [1,]  0.322 -0.407
#> [2,] -0.407  3.195

BVHAR-L

(bvhar_vhar_spec <- set_weight_bvhar(
  sigma = c(1.2, 2.3), # sigma vector
  lambda = .2, # lambda
  eps = 1e-04, # very small number
  daily = c(.5, 1), # 2-dim daily weight vector
  weekly = c(.2, .3), # 2-dim weekly weight vector
  monthly = c(.1, .1) # 2-dim monthly weight vector
))
#> Model Specification for BVHAR
#> 
#> Parameters: Coefficent matrice and Covariance matrix
#> Prior: MN_VHAR
#> ========================================================
#> 
#> Setting for 'sigma':
#> [1]  1.2  2.3
#> 
#> Setting for 'lambda':
#> [1]  0.2
#> 
#> Setting for 'eps':
#> [1]  1e-04
#> 
#> Setting for 'daily':
#> [1]  0.5  1.0
#> 
#> Setting for 'weekly':
#> [1]  0.2  0.3
#> 
#> Setting for 'monthly':
#> [1]  0.1  0.1
#> 
#> Setting for 'hierarchical':
#> [1]  FALSE

(sim_mnvhar_coef(bvhar_vhar_spec))
#> $coefficients
#>          [,1]    [,2]
#> [1,]  0.38047  0.2058
#> [2,]  0.14585  0.5599
#> [3,]  0.09613 -0.0491
#> [4,]  0.05043  0.0088
#> [5,]  0.05409  0.3732
#> [6,] -0.00569 -0.0482
#> 
#> $covmat
#>       [,1]   [,2]
#> [1,] 0.566  0.263
#> [2,] 0.263 14.995