Introduction to bvhar

library(bvhar)

Data

Looking at VAR and VHAR, you can learn how the models work and how to perform this package.

ETF Dataset

This package includes some datasets. Among them, we try CBOE ETF volatility index (etf_vix). Since this is just an example, we arbitrarily extract a small number of variables: Gold, crude oil, euro currency, and china ETF.

var_idx <- c("GVZCLS", "OVXCLS", "EVZCLS", "VXFXICLS")
etf <- 
  etf_vix %>% 
  dplyr::select(dplyr::all_of(var_idx))
etf
#> # A tibble: 905 × 4
#>    GVZCLS OVXCLS EVZCLS VXFXICLS
#>     <dbl>  <dbl>  <dbl>    <dbl>
#>  1   21.5   36.5   13.2     30.2
#>  2   21.5   35.4   12.6     28.9
#>  3   22.3   35.5   13.1     29.1
#>  4   21.6   36.6   12.8     28.5
#>  5   21.2   35.6   13.3     29.5
#>  6   21.4   34.8   13.2     29.1
#>  7   21.6   34.0   13.2     28.7
#>  8   21.1   32.6   12.8     28.0
#>  9   20.3   33.5   12.7     28.9
#> 10   19.6   33.4   12.4     28.0
#> # ℹ 895 more rows

h-step ahead forecasting

For evaluation, split the data. The last 19 observations will be test set. divide_ts() function splits the time series into train-test set.

In the other vignette, we provide how to perform out-of-sample forecasting.

h <- 19
etf_eval <- divide_ts(etf, h) # Try ?divide_ts
etf_train <- etf_eval$train # train
etf_test <- etf_eval$test # test
# dimension---------
m <- ncol(etf)

Models

VAR

This package indentifies VAR(p) model by

\[\mathbf{Y}_t = \mathbf{c}+ \boldsymbol\beta_1 \mathbf{Y}_{t - 1} + \ldots + \boldsymbol\beta_p +\mathbf{Y}_{t - p} + \boldsymbol\epsilon_t\]

where \(\boldsymbol\epsilon_t \sim N(\mathbf{0}_k, \Sigma_e)\)

var_lag <- 5

The package perform VAR(p = 5) based on

\[Y_0 = X_0 A + Z\]

where

\[ Y_0 = \begin{bmatrix} \mathbf{y}_{p + 1}^T \\ \mathbf{y}_{p + 2}^T \\ \vdots \\ \mathbf{y}_n^T \end{bmatrix}_{s \times m} \equiv Y_{p + 1} \in \mathbb{R}^{s \times m} \]

by build_y0()

and

\[ X_0 = \left[\begin{array}{c|c|c|c} \mathbf{y}_p^T & \cdots & \mathbf{y}_1^T & 1 \\ \mathbf{y}_{p + 1}^T & \cdots & \mathbf{y}_2^T & 1 \\ \vdots & \vdots & \cdots & \vdots \\ \mathbf{y}_{T - 1}^T & \cdots & \mathbf{y}_{T - p}^T & 1 \end{array}\right]_{s \times k} = \begin{bmatrix} Y_p & Y_{p - 1} & \cdots & \mathbf{1}_{T - p} \end{bmatrix} \in \mathbb{R}^{s \times k} \]

by build_design(). Coefficient matrix is the form of

\[ A = \begin{bmatrix} A_1^T \\ \vdots \\ A_p^T \\ \mathbf{c}^T \end{bmatrix} \in \mathbb{R}^{k \times m} \]

This form also corresponds to the other model. Use var_lm(y, p) to model VAR(p). You can specify type = "none" to get model without constant term.

(fit_var <- var_lm(etf_train, var_lag))
#> Call:
#> var_lm(y = etf_train, p = var_lag)
#> 
#> VAR(5) Estimation using least squares
#> ====================================================
#> 
#> LSE for A1:
#>           GVZCLS_1  OVXCLS_1  EVZCLS_1  VXFXICLS_1
#> GVZCLS     0.93290    0.0545    0.0659     -0.0346
#> OVXCLS    -0.02367    1.0047   -0.1447      0.0324
#> EVZCLS    -0.00789    0.0102    0.9810      0.0199
#> VXFXICLS  -0.03868    0.0109    0.0754      0.9328
#> 
#> 
#> LSE for A2:
#>           GVZCLS_2  OVXCLS_2  EVZCLS_2  VXFXICLS_2
#> GVZCLS     -0.0781  -0.04865    0.0829      0.0561
#> OVXCLS      0.0880   0.01207    0.2729     -0.1173
#> EVZCLS      0.0195   0.00255   -0.1071     -0.0383
#> VXFXICLS    0.0896   0.04278   -0.0691      0.0419
#> 
#> 
#> LSE for A3:
#>           GVZCLS_3  OVXCLS_3  EVZCLS_3  VXFXICLS_3
#> GVZCLS      0.0424  -0.00452  -0.03245    -0.05967
#> OVXCLS     -0.0272  -0.09144  -0.05764    -0.06255
#> EVZCLS     -0.0123   0.00864   0.08693     0.00252
#> VXFXICLS   -0.0266  -0.04810   0.00851    -0.02137
#> 
#> 
#> LSE for A4:
#>           GVZCLS_4  OVXCLS_4  EVZCLS_4  VXFXICLS_4
#> GVZCLS    -0.00793   0.01072  -0.01513      0.0616
#> OVXCLS    -0.04343  -0.00377  -0.00694      0.1445
#> EVZCLS     0.00614  -0.02278  -0.01007      0.0200
#> VXFXICLS  -0.00755  -0.05555   0.08783     -0.1025
#> 
#> 
#> LSE for A5:
#>           GVZCLS_5  OVXCLS_5  EVZCLS_5  VXFXICLS_5
#> GVZCLS      0.0728  -0.01745   -0.0886   -0.017273
#> OVXCLS      0.0104   0.07151   -0.0637    0.002018
#> EVZCLS     -0.0113   0.00581    0.0202    0.000498
#> VXFXICLS   -0.0155   0.04192   -0.0254    0.093984
#> 
#> 
#> LSE for constant:
#>   GVZCLS    OVXCLS    EVZCLS  VXFXICLS  
#>    0.571     0.145     0.129     0.875  
#> 
#> 
#> --------------------------------------------------
#> *_j of the Coefficient matrix: corresponding to the j-th VAR lag

The package provide S3 object.

# class---------------
class(fit_var)
#> [1] "varlse"   "bvharmod"
# inheritance---------
is.varlse(fit_var)
#> [1] TRUE
# names---------------
names(fit_var)
#>  [1] "coefficients"  "fitted.values" "residuals"     "covmat"       
#>  [5] "df"            "m"             "obs"           "y0"           
#>  [9] "p"             "totobs"        "process"       "type"         
#> [13] "design"        "y"             "call"

VHAR

Consider Vector HAR (VHAR) model.

\[\mathbf{Y}_t = \mathbf{c}+ \Phi^{(d)} + \mathbf{Y}_{t - 1} + \Phi^{(w)} \mathbf{Y}_{t - 1}^{(w)} + \Phi^{(m)} \mathbf{Y}_{t - 1}^{(m)} + \boldsymbol\epsilon_t\]

where \(\mathbf{Y}_t\) is daily RV and

\[\mathbf{Y}_t^{(w)} = \frac{1}{5} \left( \mathbf{Y}_t + \cdots + \mathbf{Y}_{t - 4} \right)\]

is weekly RV

and

\[\mathbf{Y}_t^{(m)} = \frac{1}{22} \left( \mathbf{Y}_t + \cdots + \mathbf{Y}_{t - 21} \right)\]

is monthly RV. This model can be expressed by

\[Y_0 = X_1 \Phi + Z\]

where

\[ \Phi = \begin{bmatrix} \Phi^{(d)T} \\ \Phi^{(w)T} \\ \Phi^{(m)T} \\ \mathbf{c}^T \end{bmatrix} \in \mathbb{R}^{(3m + 1) \times m} \]

Let \(T\) be

\[ \mathbb{C}_0 \mathpunct{:}=\begin{bmatrix} 1 & 0 & \cdots & 0 & 0 & \cdots & 0 \\ 1 / 5 & 1 / 5 & \cdots & 1 / 5 & 0 & \cdots & 0 \\ 1 / 22 & 1 / 22 & \cdots & 1 / 22 & 1 / 22 & \cdots & 1 / 22 \end{bmatrix} \otimes I_m \in \mathbb{R}^{3m \times 22m} \]

and let \(\mathbb{C}_{HAR}\) be

\[ \mathbb{C}_{HAR} \mathpunct{:}=\left[\begin{array}{c|c} T & \mathbf{0}_{3m} \\ \hline \mathbf{0}_{3m}^T & 1 \end{array}\right] \in \mathbb{R}^{(3m + 1) \times (22m + 1)} \]

Then for \(X_0\) in VAR(p),

\[ X_1 = X_0 \mathbb{C}_{HAR}^T = \begin{bmatrix} \mathbf{y}_{22}^T & \mathbf{y}_{22}^{(w)T} & \mathbf{y}_{22}^{(m)T} & 1 \\ \mathbf{y}_{23}^T & \mathbf{y}_{23}^{(w)T} & \mathbf{y}_{23}^{(m)T} & 1 \\ \vdots & \vdots & \vdots & \vdots \\ \mathbf{y}_{T - 1}^T & \mathbf{y}_{T - 1}^{(w)T} & \mathbf{y}_{T - 1}^{(m)T} & 1 \end{bmatrix} \in \mathbb{R}^{s \times (3m + 1)} \]

This package fits VHAR by scaling VAR(p) using \(\mathbb{C}_{HAR}\) (scale_har(m, week = 5, month = 22)). Use vhar_lm(y) to fit VHAR. You can specify type = "none" to get model without constant term.

(fit_har <- vhar_lm(etf_train))
#> Call:
#> vhar_lm(y = etf_train)
#> 
#> VHAR Estimation====================================================
#> 
#> LSE for day:
#>           GVZCLS_day  OVXCLS_day  EVZCLS_day  VXFXICLS_day
#> GVZCLS       0.87561      0.0447      0.1623      -0.03772
#> OVXCLS       0.04147      0.9942     -0.0605      -0.09361
#> EVZCLS       0.00305      0.0281      0.9206      -0.00748
#> VXFXICLS     0.01021      0.0569      0.0440       0.91713
#> 
#> 
#> LSE for week:
#>           GVZCLS_week  OVXCLS_week  EVZCLS_week  VXFXICLS_week
#> GVZCLS        0.01622      -0.0554      -0.1608         0.0637
#> OVXCLS       -0.07093      -0.0373       0.2000         0.1034
#> EVZCLS       -0.00334      -0.0414      -0.0101         0.0239
#> VXFXICLS     -0.03756      -0.0787      -0.0135         0.0480
#> 
#> 
#> LSE for month:
#>           GVZCLS_month  OVXCLS_month  EVZCLS_month  VXFXICLS_month
#> GVZCLS        0.084981       0.00359        0.0228         -0.0299
#> OVXCLS        0.045986       0.03825       -0.1564         -0.0157
#> EVZCLS       -0.000597       0.02030        0.0501         -0.0138
#> VXFXICLS      0.041648       0.01263        0.0639         -0.0371
#> 
#> 
#> LSE for constant:
#>   GVZCLS    OVXCLS    EVZCLS  VXFXICLS  
#>    0.491     0.135     0.105     0.926  
#> 
#> 
#> --------------------------------------------------
#> *_day, *_week, *_month of the Coefficient matrix: daily, weekly, and monthly term in the VHAR model
# class----------------
class(fit_har)
#> [1] "vharlse"  "bvharmod"
# inheritance----------
is.varlse(fit_har)
#> [1] FALSE
is.vharlse(fit_har)
#> [1] TRUE
# complements----------
names(fit_har)
#>  [1] "coefficients"  "fitted.values" "residuals"     "covmat"       
#>  [5] "df"            "m"             "obs"           "y0"           
#>  [9] "p"             "week"          "month"         "totobs"       
#> [13] "process"       "type"          "HARtrans"      "design"       
#> [17] "y"             "call"

BVAR

bvhar package provides two prior framework - Minnesota prior and flat prior.

Minnesota prior

  • Litterman (1986) and Bańbura et al. (2010)
  • All the equations are centered around the random walk with drift.
  • Prior mean: Recent lags provide more reliable information than the more distant ones.
  • Prior variance: Own lags explain more of the variation of a given variable than the lags of other variables in the equation.

First specify the prior using set_bvar(sigma, lambda, delta, eps = 1e-04).

bvar_lag <- 5
sig <- apply(etf_train, 2, sd) # sigma vector
lam <- .2 # lambda
delta <- rep(0, m) # delta vector (0 vector since RV stationary)
eps <- 1e-04 # very small number
(bvar_spec <- set_bvar(sig, lam, delta, eps))
#> 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    EVZCLS  VXFXICLS  
#>     3.77     10.63      2.27      3.81  
#> 
#> Setting for 'lambda':
#> [1]  0.2
#> 
#> Setting for 'delta':
#> [1]  0  0  0  0
#> 
#> Setting for 'eps':
#> [1]  1e-04

In turn, bvar_minnesota(y, p, bayes_spec, include_mean = TRUE) fits BVAR(p).

  • y: Multivariate time series data. It should be data frame or matrix, which means that every column is numeric. Each column indicates variable, i.e. it sould be wide format.
  • p: Order of BVAR
  • bayes_spec: Output of set_bvar()
  • include_mean = TRUE: By default, you include the constant term in the model.
(fit_bvar <- bvar_minnesota(etf_train, bvar_lag, bvar_spec))
#> Call:
#> bvar_minnesota(y = etf_train, p = bvar_lag, bayes_spec = bvar_spec)
#> 
#> BVAR(5) with Minnesota Prior
#> ====================================================
#> 
#> A ~ Matrix Normal (Mean, Precision, Scale = Sigma)
#> ====================================================
#> Matrix Normal Mean for A1 part:
#>           GVZCLS_1  OVXCLS_1  EVZCLS_1  VXFXICLS_1
#> GVZCLS      0.7771   0.00915    0.0628      0.0193
#> OVXCLS      0.0445   0.70884    0.1111      0.0167
#> EVZCLS      0.0104   0.01068    0.7036      0.0266
#> VXFXICLS    0.0214   0.00673    0.1044      0.7674
#> 
#> 
#> Matrix Normal Mean for A2 part:
#>            GVZCLS_2   OVXCLS_2  EVZCLS_2  VXFXICLS_2
#> GVZCLS     0.082726  -0.006736  -0.01387     -0.0020
#> OVXCLS     0.000568   0.140781   0.02419     -0.0436
#> EVZCLS    -0.004778   0.000769   0.11756     -0.0114
#> VXFXICLS   0.013424  -0.003779   0.00369      0.1063
#> 
#> 
#> Matrix Normal Mean for A3 part:
#>           GVZCLS_3   OVXCLS_3  EVZCLS_3  VXFXICLS_3
#> GVZCLS     0.03574  -0.003825  -0.01504    -0.00581
#> OVXCLS    -0.01733   0.054037   0.00196    -0.01874
#> EVZCLS    -0.00391   0.000223   0.05065    -0.00168
#> VXFXICLS  -0.00891  -0.006167  -0.00435     0.02176
#> 
#> 
#> Matrix Normal Mean for A4 part:
#>           GVZCLS_4   OVXCLS_4   EVZCLS_4  VXFXICLS_4
#> GVZCLS     0.02474  -0.001932  -0.011546     0.00263
#> OVXCLS    -0.00987   0.030763   0.003771     0.00845
#> EVZCLS    -0.00318  -0.000206   0.027289     0.00161
#> VXFXICLS  -0.00898  -0.004378   0.000232     0.00582
#> 
#> 
#> Matrix Normal Mean for A5 part:
#>           GVZCLS_5   OVXCLS_5  EVZCLS_5  VXFXICLS_5
#> GVZCLS     0.01986  -1.47e-03  -0.00970     0.00127
#> OVXCLS    -0.00357   2.10e-02   0.00339     0.01017
#> EVZCLS    -0.00284   1.59e-05   0.01729     0.00159
#> VXFXICLS  -0.00431  -1.70e-03   0.00263     0.01207
#> 
#> 
#> Matrix Normal Mean for constant part:
#>   GVZCLS    OVXCLS    EVZCLS  VXFXICLS  
#>   0.7271    0.3713    0.0971    1.2139  
#> 
#> 
#> dim(Matrix Normal precision matrix):
#> [1]  21  21
#> 
#> 
#> Sigma ~ Inverse-Wishart
#> ====================================================
#> IW scale matrix:
#>           GVZCLS  OVXCLS  EVZCLS  VXFXICLS
#> GVZCLS      1285     375     115       287
#> OVXCLS       375    3638     131       397
#> EVZCLS       115     131     220       126
#> VXFXICLS     287     397     126      1186
#> 
#> IW degrees of freedom:
#> [1] 887
#> 
#> 
#> --------------------------------------------------
#> *_j of the Coefficient matrix: corresponding to the j-th BVAR lag

It is bvarmn class. For Bayes computation, it also has other class such as normaliw and bvharmod.

# class---------------
class(fit_bvar)
#> [1] "bvarmn"   "normaliw" "bvharmod"
# inheritance---------
is.bvarmn(fit_bvar)
#> [1] TRUE
# names---------------
names(fit_bvar)
#>  [1] "coefficients"    "fitted.values"   "residuals"       "mn_prec"        
#>  [5] "iw_scale"        "iw_shape"        "df"              "m"              
#>  [9] "obs"             "prior_mean"      "prior_precision" "prior_scale"    
#> [13] "prior_shape"     "y0"              "design"          "p"              
#> [17] "totobs"          "type"            "y"               "prior_prec"     
#> [21] "call"            "process"         "spec"

Flat prior

Ghosh et al. (2018) provides flat prior for covariance matrix, i.e. non-informative. Use set_bvar_flat(U).

(flat_spec <- set_bvar_flat(U = 5000 * diag(m * bvar_lag + 1))) # c * I
#> Model Specification for BVAR
#> 
#> Parameters: Coefficent matrice and Covariance matrix
#> Prior: Flat
#> # Type '?bvar_flat' in the console for some help.
#> ========================================================
#> 
#> Setting for 'U':
#> # A matrix:  21 x 21 
#>        [,1]  [,2]  [,3]  [,4]  [,5]
#>  [1,]  5000     0     0     0     0
#>  [2,]     0  5000     0     0     0
#>  [3,]     0     0  5000     0     0
#>  [4,]     0     0     0  5000     0
#>  [5,]     0     0     0     0  5000
#>  [6,]     0     0     0     0     0
#>  [7,]     0     0     0     0     0
#>  [8,]     0     0     0     0     0
#>  [9,]     0     0     0     0     0
#> [10,]     0     0     0     0     0
#> # ... with 11 more rows

Then bvar_flat(y, p, bayes_spec, include_mean = TRUE):

(fit_ghosh <- bvar_flat(etf_train, bvar_lag, flat_spec))
#> Call:
#> bvar_flat(y = etf_train, p = bvar_lag, bayes_spec = flat_spec)
#> 
#> BVAR(5) with Flat Prior
#> ====================================================
#> 
#> A ~ Matrix Normal (Mean, U^{-1}, Scale 2 = Sigma)
#> ====================================================
#> Matrix Normal Mean for A1 part:
#>           GVZCLS_1  OVXCLS_1  EVZCLS_1  VXFXICLS_1
#> GVZCLS      0.3128    0.0460    0.0205      0.0457
#> OVXCLS      0.0440    0.4128    0.0186      0.0314
#> EVZCLS      0.0174    0.0230    0.1334      0.0410
#> VXFXICLS    0.0477    0.0464    0.0481      0.3197
#> 
#> 
#> Matrix Normal Mean for A2 part:
#>           GVZCLS_2  OVXCLS_2  EVZCLS_2  VXFXICLS_2
#> GVZCLS     0.19811   0.00287   0.00771      0.0194
#> OVXCLS     0.01595   0.23709   0.00978     -0.0111
#> EVZCLS     0.00425   0.01205   0.11715      0.0201
#> VXFXICLS   0.02725   0.01516   0.03513      0.2162
#> 
#> 
#> Matrix Normal Mean for A3 part:
#>           GVZCLS_3  OVXCLS_3  EVZCLS_3  VXFXICLS_3
#> GVZCLS     0.13884   -0.0153  -0.00102     0.00591
#> OVXCLS    -0.00741    0.1304   0.00361    -0.02452
#> EVZCLS    -0.00353    0.0079   0.10724     0.01150
#> VXFXICLS   0.00797   -0.0146   0.02633     0.14628
#> 
#> 
#> Matrix Normal Mean for A4 part:
#>           GVZCLS_4  OVXCLS_4  EVZCLS_4  VXFXICLS_4
#> GVZCLS     0.11392  -0.01775  -0.00653     0.00884
#> OVXCLS    -0.01626   0.09421   0.00195    -0.00572
#> EVZCLS    -0.00847   0.00565   0.10046     0.01098
#> VXFXICLS  -0.00356  -0.03067   0.02292     0.10552
#> 
#> 
#> Matrix Normal Mean for A5 part:
#>           GVZCLS_5  OVXCLS_5  EVZCLS_5  VXFXICLS_5
#> GVZCLS     0.11282   -0.0208  -0.01028     0.01136
#> OVXCLS    -0.01507    0.1004   0.00155     0.00973
#> EVZCLS    -0.01252    0.0104   0.09667     0.01361
#> VXFXICLS  -0.00492   -0.0215   0.02342     0.10353
#> 
#> 
#> Matrix Normal Mean for constant part:
#>   GVZCLS    OVXCLS    EVZCLS  VXFXICLS  
#> 5.49e-03  7.82e-04  8.83e-05  9.82e-03  
#> 
#> 
#> dim(Matrix Normal precision matrix):
#> [1]  21  21
#> 
#> 
#> Sigma ~ Inverse-Wishart
#> ====================================================
#> IW scale matrix:
#>           GVZCLS  OVXCLS  EVZCLS  VXFXICLS
#> GVZCLS      2483     595     209       535
#> OVXCLS       595    3580     216       578
#> EVZCLS       209     216     771       354
#> VXFXICLS     535     578     354      2472
#> 
#> 
#> --------------------------------------------------
#> *_j of the Coefficient matrix: corresponding to the j-th BVAR lag
# class---------------
class(fit_ghosh)
#> [1] "bvarflat" "normaliw" "bvharmod"
# inheritance---------
is.bvarflat(fit_ghosh)
#> [1] TRUE
# names---------------
names(fit_ghosh)
#>  [1] "coefficients"    "fitted.values"   "residuals"       "mn_prec"        
#>  [5] "iw_scale"        "iw_shape"        "df"              "p"              
#>  [9] "m"               "obs"             "totobs"          "process"        
#> [13] "spec"            "type"            "call"            "prior_mean"     
#> [17] "prior_precision" "y0"              "design"          "y"

BVHAR

Consider the VAR(22) form of VHAR.

\[ \begin{aligned} \mathbf{Y}_t = \mathbf{c}& + \left( \Phi^{(d)} + \frac{1}{5} \Phi^{(w)} + \frac{1}{22} \Phi^{(m)} \right) \mathbf{Y}_{t - 1} \\ & + \left( \frac{1}{5} \Phi^{(w)} + \frac{1}{22} \Phi^{(m)} \right) \mathbf{Y}_{t - 2} + \cdots \left( \frac{1}{5} \Phi^{(w)} + \frac{1}{22} \Phi^{(m)} \right) \mathbf{Y}_{t - 5} \\ & + \frac{1}{22} \Phi^{(m)} \mathbf{Y}_{t - 6} + \cdots + \frac{1}{22} \Phi^{(m)} \mathbf{Y}_{t - 22} \end{aligned} \]

What does Minnesota prior mean in VHAR model?

For more simplicity, write coefficient matrices by \(\Phi^{(1)}, \Phi^{(2)}, \Phi^{(3)}\). If we apply the prior in the same way, Minnesota moment becomes

\[ E \left[ (\Phi^{(l)})_{ij} \right] = \begin{cases} \delta_i & j = i, \; l = 1 \\ 0 & o/w \end{cases} \quad \mathrm{Var}\left[ (\Phi^{(l)})_{ij} \right] = \begin{cases} \frac{\lambda^2}{l^2} & j = i \\ \nu \frac{\lambda^2}{l^2} \frac{\sigma_i^2}{\sigma_j^2} & o/w \end{cases} \]

We call this VAR-type Minnesota prior or BVHAR-S.

BVHAR-S

set_bvhar(sigma, lambda, delta, eps = 1e-04) specifies VAR-type Minnesota prior.

(bvhar_spec_v1 <- set_bvhar(sig, lam, delta, eps))
#> 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    EVZCLS  VXFXICLS  
#>     3.77     10.63      2.27      3.81  
#> 
#> Setting for 'lambda':
#> [1]  0.2
#> 
#> Setting for 'delta':
#> [1]  0  0  0  0
#> 
#> Setting for 'eps':
#> [1]  1e-04

bvhar_minnesota(y, har = c(5, 22), bayes_spec, include_mean = TRUE) can fit BVHAR with this prior. This is the default prior setting.

(fit_bvhar_v1 <- bvhar_minnesota(etf_train, bayes_spec = bvhar_spec_v1))
#> Call:
#> bvhar_minnesota(y = etf_train, bayes_spec = bvhar_spec_v1)
#> 
#> BVHAR with Minnesota Prior
#> ====================================================
#> 
#> Phi ~ Matrix Normal (Mean, Scale 1, Scale 2 = Sigma)
#> ====================================================
#> Matrix Normal Mean for day:
#>           GVZCLS_day  OVXCLS_day  EVZCLS_day  VXFXICLS_day
#> GVZCLS        0.7808     0.00502      0.0595        0.0181
#> OVXCLS        0.0419     0.75268      0.1293       -0.0129
#> EVZCLS        0.0109     0.00957      0.7248        0.0213
#> VXFXICLS      0.0252     0.00347      0.1003        0.8042
#> 
#> 
#> Matrix Normal Mean for week:
#>           GVZCLS_week  OVXCLS_week  EVZCLS_week  VXFXICLS_week
#> GVZCLS         0.1160    -0.007669     -0.03216       -0.00181
#> OVXCLS        -0.0211     0.152822      0.02794       -0.00136
#> EVZCLS        -0.0103    -0.000372      0.13214       -0.00158
#> VXFXICLS      -0.0177    -0.010393      0.00229        0.10274
#> 
#> 
#> Matrix Normal Mean for month:
#>           GVZCLS_month  OVXCLS_month  EVZCLS_month  VXFXICLS_month
#> GVZCLS         0.05269      -0.00345       -0.0110        -0.00302
#> OVXCLS         0.00556       0.05169       -0.0225        -0.01051
#> EVZCLS        -0.00441       0.00410        0.0502        -0.00138
#> VXFXICLS       0.00733      -0.00321        0.0109         0.00573
#> 
#> 
#> Matrix Normal Mean for constant part:
#>   GVZCLS    OVXCLS    EVZCLS  VXFXICLS  
#>   0.6107    0.1410    0.0741    1.1544  
#> 
#> 
#> dim(Matrix Normal precision matrix):
#> [1]  13  13
#> 
#> 
#> Sigma ~ Inverse-Wishart
#> ====================================================
#> IW scale matrix:
#>           GVZCLS  OVXCLS  EVZCLS  VXFXICLS
#> GVZCLS      1268     366     114       286
#> OVXCLS       366    3742     131       378
#> EVZCLS       114     131     219       121
#> VXFXICLS     286     378     121      1189

This model is bvharmn class.

# class---------------
class(fit_bvhar_v1)
#> [1] "bvharmn"  "normaliw" "bvharmod"
# inheritance---------
is.bvharmn(fit_bvhar_v1)
#> [1] TRUE
# names---------------
names(fit_bvhar_v1)
#>  [1] "coefficients"    "fitted.values"   "residuals"       "mn_prec"        
#>  [5] "iw_scale"        "iw_shape"        "df"              "m"              
#>  [9] "obs"             "prior_mean"      "prior_precision" "prior_scale"    
#> [13] "prior_shape"     "y0"              "design"          "p"              
#> [17] "week"            "month"           "totobs"          "type"           
#> [21] "HARtrans"        "y"               "prior_prec"      "call"           
#> [25] "process"         "spec"

BVHAR-L

Set \(\delta_i\) for weekly and monthly coefficient matrices in above Minnesota moments:

\[ E \left[ (\Phi^{(l)})_{ij} \right] = \begin{cases} d_i & j = i, \; l = 1 \\ w_i & j = i, \; l = 2 \\ m_i & j = i, \; l = 3 \end{cases} \]

i.e. instead of one delta vector, set three vector

  • daily
  • weekly
  • monthly

This is called VHAR-type Minnesota prior or BVHAR-L.

set_weight_bvhar(sigma, lambda, eps, daily, weekly, monthly) defines BVHAR-L.

daily <- rep(.1, m)
weekly <- rep(.1, m)
monthly <- rep(.1, m)
(bvhar_spec_v2 <- set_weight_bvhar(sig, lam, eps, daily, weekly, monthly))
#> 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    EVZCLS  VXFXICLS  
#>     3.77     10.63      2.27      3.81  
#> 
#> Setting for 'lambda':
#> [1]  0.2
#> 
#> Setting for 'eps':
#> [1]  1e-04
#> 
#> Setting for 'daily':
#> [1]  0.1  0.1  0.1  0.1
#> 
#> Setting for 'weekly':
#> [1]  0.1  0.1  0.1  0.1
#> 
#> Setting for 'monthly':
#> [1]  0.1  0.1  0.1  0.1

bayes_spec option of bvhar_minnesota() gets this value, so you can use this prior intuitively.

fit_bvhar_v2 <- bvhar_minnesota(
  etf_train, 
  bayes_spec = bvhar_spec_v2
)
fit_bvhar_v2
#> Call:
#> bvhar_minnesota(y = etf_train, bayes_spec = bvhar_spec_v2)
#> 
#> BVHAR with Minnesota Prior
#> ====================================================
#> 
#> Phi ~ Matrix Normal (Mean, Scale 1, Scale 2 = Sigma)
#> ====================================================
#> Matrix Normal Mean for day:
#>           GVZCLS_day  OVXCLS_day  EVZCLS_day  VXFXICLS_day
#> GVZCLS        0.7780     0.00503      0.0605        0.0183
#> OVXCLS        0.0435     0.74834      0.1234       -0.0110
#> EVZCLS        0.0112     0.00947      0.7218        0.0211
#> VXFXICLS      0.0253     0.00363      0.0997        0.8021
#> 
#> 
#> Matrix Normal Mean for week:
#>           GVZCLS_week  OVXCLS_week  EVZCLS_week  VXFXICLS_week
#> GVZCLS         0.1180    -0.007746     -0.03245       -0.00215
#> OVXCLS        -0.0213     0.155999      0.02539       -0.00152
#> EVZCLS        -0.0104    -0.000477      0.13525       -0.00186
#> VXFXICLS      -0.0182    -0.010424      0.00136        0.10491
#> 
#> 
#> Matrix Normal Mean for month:
#>           GVZCLS_month  OVXCLS_month  EVZCLS_month  VXFXICLS_month
#> GVZCLS         0.05510      -0.00346      -0.01138        -0.00359
#> OVXCLS         0.00473       0.05518      -0.02463        -0.01064
#> EVZCLS        -0.00451       0.00393       0.05338        -0.00178
#> VXFXICLS       0.00681      -0.00319       0.00994         0.00823
#> 
#> 
#> Matrix Normal Mean for constant part:
#>   GVZCLS    OVXCLS    EVZCLS  VXFXICLS  
#>   0.5989    0.1195    0.0759    1.1240  
#> 
#> 
#> dim(Matrix Normal precision matrix):
#> [1]  13  13
#> 
#> 
#> Sigma ~ Inverse-Wishart
#> ====================================================
#> IW scale matrix:
#>           GVZCLS  OVXCLS  EVZCLS  VXFXICLS
#> GVZCLS      1252     367     114       286
#> OVXCLS       367    3607     130       380
#> EVZCLS       114     130     214       121
#> VXFXICLS     286     380     121      1174

Inference

summary() method for normaliw object performs MCMC based on conjugate MNIW distribution.

(inf_bvar <- summary(fit_bvar, num_iter = 10000, num_burn = 5000))
#> Call:
#> bvar_minnesota(y = etf_train, p = bvar_lag, bayes_spec = bvar_spec)
#> 
#> BVAR(5) with Minnesota Prior
#> ====================================================
#> Phi ~ Matrix Normal (Mean, Precision, Scale = Sigma)
#> Sigma ~ Inverse-Wishart (IW Scale, IW df)
#> 
#> 
#> Conjugate MCMC:
#> ====================================================
#> Total number of iteration: 10000
#> Number of burn-in: 5000
#> ====================================================
#> 
#> Parameter record:
#> # A draws_df: 5000 iterations, 1 chains, and 94 variables
#>     alpha[1]  alpha[2]  alpha[3]  alpha[4]  alpha[5]  alpha[6]  alpha[7]
#> 1      0.747   0.02977   0.05125  -0.01931    0.0985  -0.00565   0.04407
#> 2      0.755   0.01932   0.02100   0.04635    0.1096  -0.01795  -0.01429
#> 3      0.789   0.00994   0.09733   0.04070    0.0611  -0.00808  -0.00896
#> 4      0.764   0.01586  -0.06864  -0.00057    0.0976  -0.00903   0.02693
#> 5      0.753   0.01811   0.09808   0.01871    0.0892  -0.02147  -0.03533
#> 6      0.773   0.03325   0.00597  -0.00851    0.1269  -0.01797   0.03579
#> 7      0.825   0.01409   0.01920   0.01394    0.0619  -0.00515  -0.02075
#> 8      0.729   0.01241   0.04429   0.04073    0.0919  -0.00197  -0.07180
#> 9      0.774   0.03300  -0.03091   0.02414    0.1010  -0.00558   0.02616
#> 10     0.773  -0.00326  -0.07341   0.00432    0.0829   0.00462   0.02083
#>     alpha[8]
#> 1   -0.03839
#> 2   -0.00821
#> 3   -0.01190
#> 4   -0.00685
#> 5   -0.02286
#> 6    0.03043
#> 7   -0.02632
#> 8   -0.00278
#> 9    0.02287
#> 10   0.01300
#> # ... with 4990 more draws, and 86 more variables
#> # ... hidden reserved variables {'.chain', '.iteration', '.draw'}

autoplot() gives several Bayes visualization based on bayesplot package.

autoplot(inf_bvar, type = "trace", pars = c("alpha[1]", "alpha[2]"))