r - 事前確率全体でのモデルの限界尤度の変動が非常に小さい

Question

この記事のミネソタ州のハイパーパラメータを使用して、ベイジアン ARDL のさまざまな事前確率でモデルの限界尤度を計算しようとしています。また、正規の限界尤度の式を使用しています-Koop 2003、p.41のinvガンマモデルを使用して、異なる事前確率と生成されたデータが与えられた場合のモデルの限界尤度を計算します。コードは次のように記述します。

この部分は次のデータを生成します。

# Data
r1 <- rnorm(200)
r2 <- rnorm(200)
r3 <- rnorm(200)
e <- rnorm(200, 0, .4)

x1 <- r1 + r2
x2 <- 0.5*r1 + r3

Y <- matrix(0, 200, 1)
Y[1, 1] <- 2

for (i in 2:nrow(Y)) {
  Y[i, 1] <- 0.5 + Y[i-1, 1] * 0.7 + x1[i] * 2 + x1[i-1] * 0.5 + e[i] + 3
}
plot(Y, type = "l")
X <- as.matrix(cbind(Y[-200, ],  x1[-1], x1[-200], x2[-1], x2[-200], 1))
Y <- as.matrix(Y[-1, ])

この部分は関数を定義します。最初の関数は、 に含まれる事前に指定されたハイパーパラメーターを生成しzます。2 番目の関数は、[Koop、2003 年] からモデルの限界尤度を計算します。

ylag <-  1
xlag <-  1

  vr_mat <- matrix(0, 1, 1 + (ncol(X) - ylag - 1) / (xlag + 1))
  for (i in 1:ncol(vr_mat)) {
    if (i == 1) {
      x <- as.matrix(cbind(1, Y[-nrow(Y)]))
      vr_mat[1, i] <- sd(Y[-1, ] - x %*% solve(t(x) %*% x) %*% t(x) %*% Y[-1, ])
    } else {
      x <- as.matrix(cbind(1, X[-nrow(X), ylag + 1 + (i - 2) * (xlag + 1)]))
      y <- as.matrix(cbind(1, X[-1, ylag + 1 + (i - 2) * (xlag + 1)]))
      vr_mat[1, i] <- sd(y - x %*% solve(t(x) %*% x) %*% t(x) %*% y)
  }
}

gen_vars <- function(vr_mat, z, C) {
  x <- unlist(z)
  K  <- (ncol(X) - 1 - ylag) / (xlag + 1)
  ys <- x[1] * x[3] * (2:(1 + ylag))^-x[2]
  vars <- matrix(0, 1, K * xlag)
  L <- list()
  for (i in 1:K) {
    L[[i]] <- x[1] * (1 + 0 : xlag) ^ -x[2] * vr_mat[1, i + 1] / vr_mat[1, 1]
  }
  return(diag(c(ys, unlist(L), C)))
}

mlik <- function(beta, sig0, v0, d0) {
  s0 <- v0/d0
  bhat <- as.matrix(lm.fit(X, Y, method = "qr")$coefficients)
  v <- nrow(Y) - ncol(X)
  vbar <- v + v0
  ehat <- Y - X %*% bhat
  s <- (t(ehat) %*% ehat) / v
  vs <- (v0*s0)+(v*s) + t(bhat - beta) %*% solve(sig0 + solve(t(X) %*% X)) %*% (bhat-beta)
  c <- (gamma(vbar/2)*(v0*s0)^(v0/2)) / (gamma(v0/2)*(pi^(nrow(Y)/2)))
  VV <- solve(sig0+solve(t(X)%*%X))
  return(c * ((det(VV)/det(sig0))^0.5)*((vs)^(-vbar/2)))
}

この部分では、さまざまなハイパーパラメータで限界尤度を計算します。

v0 <- 1
d0 <- 0.1

# Check priors. Chib
beta <- c(0.7, 2, 0.5, 0, 0, 3.5)

grid <- expand.grid(z1 = seq(1, 10, by = 0.25),
                    z2 = seq(1, 10, by = 0.25),
                    z3 = seq(1, 3, by = 0.25))

res <- matrix(0, nrow(grid), 1)
for (i in 1:nrow(grid)) {
  sig0 <- gen_vars(vr_mat, grid[i, ], 5)
  res[i, 1] <- mlik(beta, sig0, v0, d0)
  cat("\n", i)
}

ただし、変動をチェックresすると、非常に小さくなり、最大 - 最小が非常に小さいことも保持されresます。このように動作するはずですか？