将方差gamma拟合到数据

Question

我使用R studio来估计Variance Gamma下数据的参数。我想将这些数据与数据进行拟合，并找到参数的估计值。我拥有的代码是

x<-c(1291,849,238,140,118,108,87,70,63,58,50,47,21,21,19)
library(VarianceGamma)
init<-c(0,0.5,0,0.5)
vgFit(x, freq = NULL, breaks = NULL, paramStart = init, startMethod = "Nelder-Mead", startValues = "SL", method = "Nelder-Mead", hessian = FALSE, plots = TRUE)

我得到的错误是：

optim错误(paramStart，llsklp，NULL，method = startMethodSL，hessian = FALSE，：函数不能在初始参数求值我不确定问题是什么？

Answer 1

该错误可能表明存在分歧。根据你之前的问题，我猜x是股票价值的原始数字。因此，在对每个时间单位的变化进行建模之前，可能需要一个log-transformation (例如，每日回报)。

x <- c(1291,849,238,140,118,108,87,70,63,58,50,47,21,21,19)
dx <- log(x)[2:length(dx)] - log(x)[1:(length(dx)-1)]

vgFit(dx)
#Parameter estimates:
#     vgC     sigma     theta        nu  
# 0.16887   0.03128  -0.47164   0.27558

我们可能想要与模拟数据进行比较。我实现了两个方法，对于大量的观察值nt，它们似乎是等效的。

方法二参照如下：

​

#Simulating VG as a time-fixed Brownian Motion
set.seed(1) 

nt = 15 #number of observations
T = nt - 1 #total time
dt = rep(T/(nt-1), nt-1) #fixed time increments
r = 1 + 0.16887 #interest rate
vgC = (r-1)
sigma = 0.03128 
theta = -0.47164 
nu = 0.27558

V_ = rep(NA,nt) #Simulations for log stock value
V_[1] = 7.163172 #log(x[1])
V2_ = V_ #alternative simulation method
for(i in 2:nt) 
    {#method 1: by VarianceGamma package
    V_[i] <- V_[i-1] + rvg(1,vgC=vgC*dt[i-1], sigma=sigma, theta=theta, nu=nu)
    
    #method 2: by R built-in packages
    gamma_i<-rgamma(1, shape=dt[i-1]/nu, scale = nu)
    normal<-rnorm(1, mean=0, sd=sigma*sqrt(gamma_i))
    V2_[i] <- V2_[i-1] + vgC*dt[i-1] + theta*gamma_i + normal
    }

# Visual comparison
x11(width=4,height=4)
plot(x, xlab='Time',ylab='Stock value',type='l')
lines(exp(V_), col='red')
lines(exp(V2_), col='blue')
legend('topright',legend=c('Observed','Method1','Method2'),fill=c('black','red','blue'))

由此得到的参数表明，由于样本大小较小的nt，估计值不稳定

#The real parameter:
c(vgC*dt[1], sigma, theta, nu).
#     vgC     sigma     theta        nu 
# 0.16887   0.03128  -0.47164   0.27558

#Parameter estimates for 1st data set:
dV = V_[2:nt] - V_[1:(nt-1)]
vgFit(dV)
#    vgC    sigma    theta       nu  
#-0.9851   0.3480   1.2382   2.0000

#Parameter estimates for 2nd data set:
dV2 = V2_[2:nt] - V2_[1:(nt-1)]
vgFit(dV2)
#     vgC     sigma     theta        nu  
#-0.78033   0.07641   0.52414   0.11840

此外，rvg函数假定时间增量是固定的。我们可以使用log-likelihood方法通过@Louis Marascio的answer来放松这一假设。

#Simulating VG as a time-changed Brownian Motion
set.seed(1) 

nt = 100 #Increase the number of observations!
T = nt-1
dt = runif(nt-1) #random time increments 
dt = dt/sum(dt)*T
r = 1 + 0.16887
vgC = (r-1)
sigma = 0.03128 
theta = -0.47164 
nu = 0.27558

V_ = rep(NA,nt) #simulations for log stock value
V_[1] = 7.163172
for(i in 2:nt) 
    {V_[i] <- V_[i-1] + rvg(1,vgC=vgC*dt[i-1], sigma=sigma, theta=theta, nu=nu)
    }
dV = V_[2:nt] - V_[1:(nt-1)]

# -log-likelihood function with different time increments
ll = function(par){
if(par[2]>0 & par[4]>0)
    {tem = 0
    for (i in 1:(length(dV)))
        {tem = tem - log(dvg(dV[i], vgC = par[1]*dt[i], sigma=par[2], theta=par[3], nu = par[4]))
        }
    return (tem)
    }
else return(Inf)}

实际上，通过放松固定时间假设，结果显示了更好的估计：

#The real parameters:
c(vgC, sigma, theta, nu)
#       vgC      sigma      theta         nu
#   0.16887    0.03128   -0.47164    0.27558

#Assuming fixed time increments
vgFit(dV)$param*c(1/mean(dt),1,1,1)
#       vgC      sigma      theta         nu 
#-0.2445969  0.3299023 -0.0696895  1.5623556

#Assuming different time increments
optim(vgFit(dV)$param*c(1/mean(dt),1,1,1),ll,
    method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Brent")[5])
#        vgC       sigma       theta          nu 
# 0.16503125  0.03241617 -0.50193694  0.28221985

最后，可以通过多次模拟来获得估计参数的置信区间：

set.seed(1)
out = NULL
for (j in 1:100) #100 simulations
    {V_ = rep(NA,nt)
    V_[1] = 7.163172
    for(i in 2:nt) 
        {V_[i] <- V_[i-1] + rvg(1,vgC=vgC*dt[i-1], sigma=sigma, theta=theta, nu=nu)
        }
    dV = V_[2:nt] - V_[1:(nt-1)]
    
    #to skip divergence
    tem <- try(vgFit(dV)$param)
    if (inherits(tem, "try-error")) next

    out = rbind(out,tem)
    }

apply(out,2,mean)
#       vgC      sigma      theta         nu 
#-0.8735168  0.1652970  0.4737270  0.9821458
apply(out,2,sd)
#      vgC     sigma     theta        nu 
#2.8935938 0.3092993 2.6833866 1.3161695