Simulating Meixner using Johnson translation system

[raphill]

Hello,

i want to draw random samples from the Meixner distribution using Matlab by utilizing the Johnson system for unbounded distributions. Therefore i need to estimate the parameters of the johnson system and consequently use rejection sampling to draw random samples from the Meixner distribution. However using moment estimators i end up with a mess of parameters including mainly only imaginary parts. I base my estimation on the following article by Matteo Grigoletto and Corrado Provasi: "Simulation and Estimation of the Meixner Distribution". The moments (m1 to m4) are -2.0821e-04, 8.5777e-04, -2.2295e-06 and 5.2977e-06 respectively. Could someone help me on this topic?

Thank you for your help

 

[Guaramy]

Hi raphill,

I have R-Code to simulate the Meixner Process that i can provide to you. The only problem is that i dont use the Johnson method to draw random samples for the Meixner distribuition, my plots are based in the density function.

Best Regards

 

[Shantanu Kumar]

Problem could be with moment based estimation.They are not very reliable in case of smaller samples.If possible try using MLE based estimation.

 

[raphill]

Thank you for the responses. I was able to implement an algorithm to draw random samples for the Meixner distribution. Guaramy could you provide me with your code so i can verify if your code leads to comparable results?

 

[gigisplillo]

HI Guarami. How are you? I am a writing thesis on Lèvy processes, in particular i'm constructing an hedging portfolio using the Variance Gamma process, the CGMY process and finally the Meixner process with montecarlo simulation and the FFT.
Could you help me in wrinting the R code of the Meinxer process?

 

[Guaramy]

Hi gigisplillo. I have R code that simulates and calculates a price of an European Option. I dont use monte carlo simulation, my option price equation was obtain via the neutral risk maeasure. I can send you the code without any problem, just tell me if you are interested.

 

[gigisplillo]

Of course i'm interested in.

 

[gigisplillo]

Guaramy it is possible with your R code use the fast fourier transformation?

 

[Guaramy]

I will provide you the code today, sorry for the delay but i am not living in my home country. You want to build the option price equation using the FFT? Yes, but you will have to adapt the code. Just take the part where i build the pdf of the Meixner distribution and build your formula using FFT. I hope to post the code ASAP. Sorry about the delay.

 

[Guaramy]

gigisplillo, this is my code to generate the pdf and calculate an option price via neutral risk measure. In order to calculate the complex gamma function you have to install a package which i am almost sure that is called asianoptions but confirm first please. In the first function the pdf is generated. a,b,m and d you have to estimate using the temporal series of the underlying. You can use whatever method you prefer to estimate the parameters. I personally have used in my thesis method of moments. Hope that this help you, any more question about models with using levy process don't hesitate to ask, i am not an expert but i have made a master thesis on it.

Meixner.pdf = function(x,a,b,m,d)
{

A = (2*cos(b/2))^(2*d)
B = 2*a*pi*gamma(2*d)
C = (b*(x-m))/a
D = d +(complex(1,0,1)*(x-m))/a
E = abs(cgamma(D))^2

f = (A/B)*exp(C)*E

#plot(x,f,type = "l",col="blue",lwd=2,main = "Função densidade de probabilidade de Meixner",xlab="x",ylab="f(x)")

}

Call.M = function(a,b,d,m,K,S0,theta,T)
{
c = log(K/S0)

tempG.fn=function(x)
{
b1 = a*(theta+1)+b
d1
G = Meixner.pdf (x,a,b1,d*T,m*T)

G
}

tempI.fn=function(x)
{
I = Meixner.pdf (x,a,a*theta+b,d*T,m*T)

I
}


intG.fn = function(t){sapply(t,FUN=tempG.fn)}
intI.fn = function(t){sapply(t,FUN=tempI.fn)}

IG = integrate(intG.fn,lower=c,upper=Inf,rel.tol=1e-10,subdivisions=1000000)
tempG = ifelse(IG$message == "OK",IG$value,NA)
II = integrate(intI.fn,lower=c,upper=Inf,rel.tol=1e-10,subdivisions=1000000)
tempH = ifelse(II$message == "OK",II$value,NA)

}

Call.M(a,b,d,m,K,S0,theta,T)

 

[gigisplillo]

Thank you very much Guarami. Don't worry for the delay in your reply.

 

[gigisplillo]

Good morning. How are you? Would you mind to tell based on what research paper(s) have you created the script of the Meixner process that you kindly sent me please?

 

[gigisplillo]

Good Afternoon Guarami. Could you tell me how to calculate the Meixner parameters using the method of moments through R please? Thank you in advance

 

[Guaramy]

Hey gigisplillo, I did all the code and calculation based on the
wim Schoutens article, here is the link https://perswww.kuleuven.be/~u0009713/004wsreport.pdf .
The equations to estimate the parameters are on the article

 

[gigisplillo]

Hi Guarami. I have this article. But can't write the code. I'm not an R expert. Is it possible for you to help me? I have constracted the code but I have problem with the parameters.
q<-2000 #how many variables

a<-.25 #parameter a

b<-.002 #parameter b

delta<-2 #parameter delta

ep<-0.001

A<-b/a

C<-pi/a

U<-runif(q)

Y<-ep/U

W<-runif(q)

z<-a*delta*sqrt(2*ep/pi)

g<-NULL

for(n in -10:10){for(i in 1:length(Y)){

g<-sum((-1)ˆn*exp(-nˆ2*piˆ2/(2*Cˆ2ˆY)))*exp(-Aˆ2*Y/2)}

v<-NULL

for(j in 1:length(g)){

if (g[j]>W[j]) v[j]=1

else v[j]=0}

S<-sum(Y*v)

tau<-z+S

X<-A*tau+rnorm(q)*sqrt(tau)

T<-10

#s<-0.001

t<-seq(from=0, to=T, by=T/(q-1)) # time-grid

l<-length(t)

M<-NULL

for (m in 1:l){

M[m]<-cumsum(X)[m]}

plot(t, M, type=’p’,cex=0.35,pch=18,xlab="time",ylab="M(t)")

 

[gigisplillo]

Hi Guarami sorry if I am disturbing you.
I have implemented the code for pricing a European Call Option through the Meixner process using the Montecarlo approach (this works) and the Fast Fourier Transform ( does not work).
But the problem is only to estimate the parameters of the process.
Could you give me a hand please?
I read the paper you sent me a few days ago. But I do not know to write the code to implement the method of moments.
To be fair I send you my code. This code is based on the paper attached pag 41-44.

See you soon

THE CODE:
############################Madan and Yor Simulation###############################

q <- 2000#how many variables

a <- 4.125#parameter a

b <- .0002#parameter b

delta <- 0.2#parameter delta

ep <- 0.001

A <- b/a

C <- pi/a

U <- runif(q)

Y <- ep/U

W <- runif(q)

z <- a*delta*sqrt(2*ep/pi)

g <- NULL

for(n in -10:10){

for(i in 1:length(Y)){

g <- sum((-1)^n*exp(-n^2*pi^2/(2*C^2^Y)))*exp(-A^2*Y/2)

}

}

v <- NULL

for(j in 1:length(g)){

if (g[j]>W[j]) {

v[j]=1

} else {

v[j]=0

}

}

S <- sum(Y*v)

tau <- z+S

X <- A*tau+rnorm(q)*sqrt(tau)

T <- 36

#s <- 0.001

t <- seq(from=0, to=T, by=T/(q-1)) # time-grid

l <- length(t)

M <- NULL

for (m in 1:l){

M[m] <- cumsum(X)[m]

}

plot(t, M, type='p',cex=0.35,pch=18,xlab="time",ylab="M(t)")



##################################Parameters#########################################

S0 <- 45.56

K <- 37

T <- 36

r <- 0.11

######################## Montecarlo con M##################################

require(Runuran)

distr <- udmeixner(alpha=4.125, beta=0.002, delta=2, mu=0)

gen <- pinvd.new(distr)

R <- ur(gen,n)

N <- rnorm(n,0,1)

X <- A*tau+rnorm(q)*sqrt(tau)

S <- S0 * exp (r * T + X)

payoff <- sapply(S, function (x)max(x-K,0))

mean(payoff)*exp(-r*T)

#########################FFT method con M###########################

FFTcall.price <- function(phi, S0, K, r, T, alpha = 1, N = 2^12, eta = 0.25) {

m <- r - log(phi(-(0+1i)))

phi.tilde <- function(u) (phi(u) * exp((0+1i) * u * m))^T

psi <- function(v) exp(-r * T) * phi.tilde((v - (alpha + 1) * (0+1i)))/(alpha^2 + alpha - v^2 + (0+1i) * (2 * alpha + 1) * v)

lambda <- (2 * pi)/(N * eta)

b <- 1/2 * N * lambda

ku <- -b + lambda * (0N - 1))

v <- eta * (0N - 1))

tmp <- exp((0+1i) * b * v) * psi(v) * eta * (3 + (-1)^(1:N) - ((1:N) - 1 == 0))/3

ft <- fft(tmp)

res <- exp(-alpha * ku) * ft/pi

inter <- spline(ku, Re(res), xout = log(K/S0))

return(inter$y * S0)

}

phimeixner <- function(u) {

tmp <- cos(b/2)/(cosh(a*u-b*(0+1i))/2)^(2*delta)

tmp <- tmp * exp((0+1i) * u * log(S0) + u*(0+1i)) }

FFTcall.price(phimeixner, S0 = S0, K = K, r = r, T = T)

 

[gigisplillo]

I attach the code. So you can read it better. The function FFTcall.price does not work given the random parameters.

 

[MaggieMegs]

Hi everyone!
I have the same problem: I am not able to write the R code to estimate the parameters of the Meixner Process. Does anyone have an R code to do that?
Thanks in advance!!!