0

我正在尝试使用本文中的超参数 a-la Minnesota 计算模型在贝叶斯 ARDL 中不同先验的边际可能性。我还使用来自 Koop 2003, p.41 的正常边际似然公式 - inv gamma model 来计算我的模型在不同先验和生成数据的情况下的边际似然。代码编写如下：

这部分生成数据：

# 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。第二个函数计算来自 [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会非常小，并且它还认为 max - min inres非常小。它应该像这样工作吗？

4

0 回答 0