Мостовая выборка методом Монте-Карло в 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 для оценки параметров.

Обратите внимание, что временной индекс начинается с 1 из-за синтаксиса R. Для стандартного отклонения в 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

Еще раз, Райан, я нашел другой алгоритм для выборки мостов, который попробовал самостоятельно, но я не уверен, что мои ответы верны. Я добавил свой код, вывод и алгоритм ниже, а также вывод, который, как мне кажется, должен выглядеть? Я использовал формат, аналогичный вашему коду:

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 as Variance gamma is a jump process with finite activity, the path should look like this: , 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: введите описание изображения здесь