compute_variance_decompositions R vs cpp

[adamwang15]

Hey @donotdespair

compute_variance_decompositions gives a different result from the R function compute_fevd

Here is a reproducible example

library(bsvarSIGNs)

# R code for fevd
compute_fevd = function(irf, horizon = 40) {
  N    = dim(irf)[1]
  S    = dim(irf)[4]
  fevd = array(NA,c(N,N,horizon,S))
  
  for (s in 1:S){
    
    for (i in 1:(horizon)){
      
      for (n in 1:N){
        for (nn in 1:N){
          fevd[n,nn,i,s]  = sum(irf[n,nn,1:i,s]^2)
        }
      }
      fevd[,,i,s]         = diag(1/apply(fevd[,,i,s],1,sum))%*%fevd[,,i,s]
    }
  }
  
  fevd        = 100*fevd
  fevd
}

data(optimism)

zero_irf          = matrix(0, nrow = 5, ncol = 5)
zero_irf[1, 1]    = 1
sign_irf          = array(0, dim = c(5, 5, 1))
sign_irf[2, 1, 1] = 1

specification = specify_bsvarSIGN$new(
  optimism*100,
  p        = 4,
  sign_irf = sign_irf,
  zero_irf = zero_irf
)

posterior = estimate(specification, S = 100)
irf       = compute_impulse_responses(posterior, horizon = 40)

fevd_R    = compute_fevd(irf, horizon = 10)
fevd_cpp  = compute_variance_decompositions(posterior, horizon = 10)

fevd_R[,,1,1]
fevd_cpp[,,1,1]

Seems like the R function is correct, have a look! Thanks.

Best,
Adam

[donotdespair]

Hey @adamwang15

Thank you so much for this feedback! I'm on it. This issue is given a high priority and I will address it before submitting the new package version to CRAN.

Gretings,

Tomasz

[donotdespair]

todos

• create a cpp function forecast_sigma2_sv
• create a cpp function forecast_sigma2_msh
• update the cpp forecast_bsvar_sv functions
• rewrite bsvars_fevd function to bsvars_fevd_homosk
• correct function bsvars_fevd_homosk to accommodate feedback from @adamwang15
• create a cpp function bsvars_fevd_heterosk
• update R methods for compute_variance_decompositions
• the compute_variance_decompositions methods for heteroskedastic models should be computed for the last point in the sample
• documentation
• tests
• NEWS

[donotdespair]

Hey @adamwang15

Can you have a look? It seems to me both functions, mine and yours, are fine. When I run the following script, all the values at the outcome are the same. Run the following script from a .cpp file in RStudio:

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp;
using namespace arma;


// [[Rcpp::export]]
cube bsvars_fevd_homosk (
    cube&    posterior_irf   // output of bsvars_irf
) {
  
  const int       N = posterior_irf.n_rows;
  const int       horizon = posterior_irf.n_slices;
  
  cube            aux_fevds(N, N, horizon);  // + 0 and inf horizons
  
    for (int h=0; h<horizon; h++) {
      for (int n=0; n<N; n++) {
        for (int nn=0; nn<N; nn++) {
          aux_fevds.subcube(n, nn, h, n, nn, h) = accu(square(posterior_irf.subcube(n, nn, 0, n, nn, h)));
        }
      }
      aux_fevds.slice(h)  = diagmat(1/sum(aux_fevds.slice(h), 1)) * aux_fevds.slice(h);
    }
    aux_fevds            *= 100;
    // fevds(s)              = aux_fevds;
  
  return aux_fevds;
} // END bsvars_fevd_homosk

// [[Rcpp::export]]
arma::field<arma::cube> bsvars_fevd_homoskk (
    arma::field<arma::cube>&    posterior_irf   // output of bsvars_irf
) {
  
  const int       N = posterior_irf(0).n_rows;
  const int       S = posterior_irf.n_rows;
  const int       horizon = posterior_irf(0).n_slices;
  
  field<cube>     fevds(S);
  cube            aux_fevds(N, N, horizon);  // + 0 and inf horizons
  
  for (int s=0; s<S; s++) {
    for (int h=0; h<horizon; h++) {
      for (int n=0; n<N; n++) {
        for (int nn=0; nn<N; nn++) {
          aux_fevds.subcube(n, nn, h, n, nn, h) = accu(square(posterior_irf(s).subcube(n, nn, 0, n, nn, h)));
        }
      }
      aux_fevds.slice(h)  = diagmat(1/sum(aux_fevds.slice(h), 1)) * aux_fevds.slice(h);
    }
    aux_fevds            *= 100;
    fevds(s)              = aux_fevds;
  } // END s loop
  
  return fevds;
} // END bsvars_fevd_homosk

// [[Rcpp::export]]
field<cube> a2q (
  cube& a,
  cube& b
) {
  field<cube> q(2);
  q(0) = a;
  q(1) = b;
  return q;
}

/*** R

compute_fevd = function(irf) {
  N    = dim(irf)[1]
  horizon   = dim(irf)[2]
  fevd = array(NA,c(N,N,horizon))
  
  # for (s in 1:S){
    
    for (i in 1:(horizon)){
      
      for (n in 1:N){
        for (nn in 1:N){
          fevd[n,nn,i]  = sum(irf[n,nn,1:i]^2)
        }
      }
      fevd[,,i]         = diag(1/apply(fevd[,,i],1,sum))%*%fevd[,,i]
    }
  # }
  
  fevd        = 100*fevd
  fevd
}

compute_fevds = function(irf, horizon = 2) {
  N    = dim(irf)[1]
  S    = dim(irf)[4]
  fevd = array(NA,c(N,N,horizon,S))
  
  for (s in 1:S){
    
    for (i in 1:(horizon)){
      
      for (n in 1:N){
        for (nn in 1:N){
          fevd[n,nn,i,s]  = sum(irf[n,nn,1:i,s]^2)
        }
      }
      fevd[,,i,s]         = diag(1/apply(fevd[,,i,s],1,sum))%*%fevd[,,i,s]
    }
  }
  
  fevd        = 100*fevd
  fevd
}

aa = array(1:8, c(2,2,2))
aa
bsvars_fevd_homosk(aa)
aaaa = a2q(aa,aa)
qq=bsvars_fevd_homoskk(aaaa)
qq[1]
qq[2]
compute_fevd(aa)
aaa = array(NA,c(dim(aa),2))
aaa[,,,1] = aa
aaa[,,,2] = aa
compute_fevds(aaa)
*/

What do you think?

[adamwang15]

Hey @donotdespair

Great! The two functions give identical results on my laptop as well, thanks for fixing this!