R studio中用于方差伽玛的桥采样蒙特卡罗方法

Question

我正在尝试在 R Studio 中使用桥接采样来模拟方差伽马过程的路径。我的代码是：

sigma = 0.5054
theta = 0.2464 
nu = 0.1184 
mu=1
N=2^(k)
k=5
V_<-rep(NA,252)
V_[0]<-0
G_[N]<-rgamma(1, shape=N*1/nu, scale=nu)
G_<-0
V<-rnorm(theta*G[N],sigma^2*G[N])
for(l in 1:k){
n<-2^(k-l)
for(j in 1:2^i-1){
i<-(2*j-1)*n
d1<-(n)*mu^2/nu
d2<-(n)*mu^2/nu
Y<-rbeta(1,d1,d2)
G_[i]<-G_[i-1]+(G[i+n]-G[i-n])*Y
G[i]
print(G_[i])
Z<-rnorm(0,(G_[i+n]-G_[i])*sigma^2*Y)
V_[i]<-Y*V_[i+n]+(1-Y)*V_[i-n]+Z
print(V_[i])
}
}
ts.plot(V[i])

我不确定我做错了什么。我试图遵循的算法如下图所示：

Answer 1

根据您的代码，模拟了一个数字序列。并且可以通过VarianceGamma::vgFit估计参数进行粗略的验证。

请注意，由于 R 语法，时间索引从 1 开始。方差 sqrt 用于rnorm 中的标准差。最后我可能不应该添加由于利率vgC 而导致的变化，因为它不包含在您的算法中。如果没有意义，请设置为0。

布朗桥模拟：

# Brownian-Gamma Bridge Sampling (BGBS) of a VG process
set.seed(1) 
M <- 10
nt <- 2^M + 1 #number of observations
T <- nt - 1 #total time
T_ <- seq(0, T, length.out=nt) #fixed time increments

#random time increments
#T_ = c(0, runif(nt-2), 1)
#T_ = sort(T_) * T

r <- 1 + 0.2 #interest rate
vgC <- (r-1)
sigma <- 0.5054
theta <- 0.2464 
nu <- 0.1184

V_ <- G_ <- rep(NA,nt)
V_[1] <- 0
G_[1] <- 0
G_[nt] <- rgamma(1, shape=T/nu, scale=nu)
V_[nt] <- rnorm(1, theta*G_[nt], sqrt(sigma^2*G_[nt]))

for (k in 1:M)
  {
  n <- 2^(M-k)
  for (j in 1:2^(k-1))
    {
    i <- (2*j-1) * n
    Y <- rbeta(1, (T_[i+1]-T_[i-n+1])/nu, (T_[i+n+1]-T_[i+1])/nu)
    G_[i+1] <- G_[i-n+1] + (G_[i+n+1] - G_[i-n+1]) * Y
    Z <- rnorm(1, sd=sqrt((G_[i+n+1] - G_[i+1]) * sigma^2 * Y))
    V_[i+1] <- Y * V_[i+n+1] + (1-Y) * V_[i-n+1] + Z
    }
  }
V_ <- V_ + vgC*T_ # changes due to interest rate

plot(T_, V_)

结果与估计大致相符：

#Estimated parameters:
library(VarianceGamma)
dV <- V_[2:nt] - V_[1:(nt-1)]
vgFit(dV)
>    vgC   sigma   theta      nu  
> 0.2996  0.5241  0.1663  0.1184

#Real parameters:
c(vgC, sigma, theta, nu)
>    vgC   sigma   theta      nu  
> 0.2000  0.5054  0.2464  0.1184

编辑

正如您所评论的，还有另一种类似的算法，可以以类似的方式实现。

您的代码可以修改如下：

set.seed(1) 
M <- 7
nt <- 2^M + 1
T <- nt - 1
T_ <- seq(0, T, length.out=nt)
sigma=0.008835
theta= -0.003856 
nu=0.263743  
vgc=0.004132

V_ <- G_ <- rep(1,nt)
G_[T+1] <- rgamma(1, shape=T/nu, scale=nu) #
V_[T+1] <- rnorm(1, theta*G_[T+1], sqrt(sigma^2*G_[T+1])) #
V_[1] <- 0
G_[1] <- 0
for (m in 1:M){ #
Y <- rbeta(1,T/(2^m*nu), T/(2^m*nu))
for (j in 1:2^(m-1)){ #
i <- (2*j-1)
G_[i*T/(2^m)+1] = G_[(i-1)*T/(2^m)+1]+(-G_[(i-1)*T/(2^m)+1]+G_[(i+1)*T/(2^m)+1])*Y #
b=G_[T*(i+1)/2^m+1] - G_[T*(i)/2^m+1] #
Z_i <- rnorm(1, sd=b*sigma^2*Y)
#V_[i] <- Y* V_[i+1] + (1-Y)*V_[i-1] + Z_i
V_[i*T/(2^m)+1] <- Y* V_[(i+1)*T/(2^m)+1] + (1-Y)*V_[(i-1)*T/(2^m)+1] + Z_i
 } 
 }
 V_ <- V_ + vgc*T_
 V_
 ts.plot(V_, main="BRIDGE", xlab="Time increment")

Answer 2

Ryan 再次，我发现了另一种我自己尝试过的桥接采样算法，但我不相信我的答案是正确的。我在下面添加了我的代码、输出和算法，以及我认为它应该看起来像的输出？我使用了与您的代码类似的格式：

set.seed(1) 
M <- 7
nt <- 2^M + 1 #number of observations
T <- nt - 1 #total time
T_ <- seq(0, T, length.out=nt) #fixed time increments
sigma=0.008835
theta= -0.003856 
nu=0.263743  
vgc=0.004132  
V_ <- G_ <- rep(1,nt)
G_[T] <- rgamma(1, shape=T/nu, scale=nu)
V_[T] <- rnorm(1, theta*G_[T], sqrt(sigma^2*G_[T]))
V_[1] <- 0
G_[1] <- 0
for (m in 2:M){
Y <- rbeta(1,T/(2^m*nu), T/(2^m*nu))
for (j in 2:2^(m-1)){
i <- (2*j-1)
G_[i*T/(2^m)] = G_[(i-1)*T/(2^m)]+(G_[(i-1)*T/(2^m)]+G_[(i+1)*T/(2^m)])*Y
b=G_[T*(i)/2^m] - G_[T*(i-1)/2^m]
Z_i <- rnorm(1, sd=b*sigma^2*Y)
V_[i] <- Y* V_[i+1] + (1-Y)*V_[i-1] + Z_i
 } 
 }
 V_ <- V_ + vgc*T_ # changes due to interest rate
 V_
 ts.plot(V_, main="BRIDGE", xlab="Time increment")

但是，这就是我从输出中得出的情节，如图 1 所示：

由于方差 gamma 是一个具有有限活动的跳跃过程，路径应如下所示：，这只是来自谷歌的方差 gamma 路径图像，顺序采样看起来像这样，我的目标是将其与模拟路径的桥采样进行比较。但我的输出看起来真的很不一样。请让我知道你的想法。如果我的代码有问题，请告诉我，谢谢。这是它的算法，与上面的非常相似，但略有不同：