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 开始。方差的平方根用于 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 所示：

Bu因为Variance gamma是一个有限的跳跃过程activity，路径应该是这样的：, this is just an image from google for variance gamma paths, the sequential sampling one looks like this and my aim is to compare it to Bridge sampling for simulating paths. But my output looks really different. Please let me know your thoughts. If there is an issue in my code let me know thanks. Here is algortihm for it, much similar to the one above but slightly different:

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

Bridge sampling Monte-carlo method in R studio for variance gamma

r

montecarlo