Fit Variance gamma to data

Question

I am using R studio to estimate paramters for data under Variance Gamma. I want to fit this data to the data and find estimates of parameters. The code I have is

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)

The error I got was:

Error in optim(paramStart, llsklp, NULL, method = startMethodSL, hessian = FALSE, : function cannot be evaluated at initial parameters I am not sure what the issue is?

Thank you did your reply but im still confused , my X is already a data vector ? If it's possible can you re-write it for me ? It will be much appreciated. Thank you — user13815972

Ryan SY Kwan Ryan SY Kwan · Accepted Answer · 2020-06-30T03:20:39

The error might suggest divergence. Based on your previous questions, I'm wildly guessing the x is the raw number of the stock values. So a log-transformation may be necessary before modelling the change per time unit (ex. daily returns).

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

We may want to compare with simulated data. I implemented two methods and they seem equivalent for large observation number nt.

Method 2 is according to below:

#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'))

The resulting parameters suggest unstable estimations due to small sample size 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

In addition, the rvg function is assuming fixed time increments. We can relax that hypothesis by @Louis Marascio's answer using log-likelihood approach.

#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)}

Indeed, the results show better estimation by relaxing the fixed time assumption:

#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

Finally, the confidence intervals for the estimated parameters may be obtained by multiple simulations:

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

Fit Variance gamma to data

1 Answers