R studio 中桥采样蒙特卡罗方法的方差伽马 [英] Bridge sampling Monte-carlo method in R studio for variance gamma

问题描述

我正在尝试在 R studio 中使用桥采样来模拟方差伽玛过程的路径.我的代码是:

I am using trying to use bridge sampling in R studio to simulate paths for the variance gamma process. My code is:

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])

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

I'm not sure what I've done wrong. The algorithm I am trying to follow is as below in the picture:

推荐答案

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

Based on your code, a numerical sequence was simulated. And it can be roughly validated by using VarianceGamma::vgFit to estimate the parameters.

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

Note that the time index starts from 1 due to R syntax. The sqrt of variance was used for the standard deviation in rnorm. And I probably shouldn't add the change due to interest rate vgC in the end, since it is not included in your algorithm. Please set it as 0 if it doesn't make sense.

布朗桥模拟:

# 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

编辑

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

As you commented, there is another similar algorithm and can be implemented in a similar way.

您的代码可以修改如下:

Your code could be modified as below:

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")

这篇关于R studio 中桥采样蒙特卡罗方法的方差伽马的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！