SERVING THE QUANTITATIVE FINANCE COMMUNITY

• 1
• 2

czar3k
Topic Author
Posts: 11
Joined: October 9th, 2017, 12:34 pm

### Variance Gamma Scaled Self Decomposable process (VGSSD model)

I’m using a Variance Gamma model to price some options with Monte Carlo method.
My model is calibrated to daily returns data.
There is no issue with this model as long as options are not path-dependent.
For european options  I can simulate prices of a stock only at maturity (say 90 days) by simply scaling sigma by the square root of the time to maturity  and drift by time to maturity. Thus the distribution of 90d returns remains skewed and fat tailed).
Unfortunately when daily returns need to be simulated in order to keep track of the prices’ paths the cumulated daily returns (final maturity returns) become Gaussian (loose their kurtosis and skewness).
As far as I understand from the paper by prof. C. O'Sullivan The Variance Gamma Scaled Self-Decomposable Process in Actuarial Modelling” VGSSD model overcomes these problems.
I assume that it allows for simulation of daily returns paths which when summed up to long-term returns retain kurtosis and skewness of the daily returns.
I’ve searched the web for the procedure of simulation for VGSSD model returns but in vain (only VG model simulation is available).
Below is a piece of Matlab code used for VG paths generation.
I would be grateful if you could suggest a change that should be made to convert the code to obtain VGSSD paths.
Let the extra  parameter (gamma) value be  0.5.

N_Sim = 365; // no of days to maturity
T = 1000 //no of paths
dt = 1; // step (1 = 1d)
params = [0.000164383561644   0.010468478451804  -0.000082191780822   0.800000000000000]
S0 = 1

function[S] = VGSimulate(N_Sim, T, dt, params, S0)
theta = params(3);
nu = params(4);
sigma = params(2);
mu = params(1);
g = gamrnd (1/nu ,nu ,N_Sim ,T);
randomnums = randn(N_Sim ,T);
S = S0* cumprod ([ ones(1,T); exp((mu + theta)*dt + theta*g +  sigma*sqrt(g*dt).*randomnums + dt/nu*log(1 - nu*theta/dt -  nu*sigma^2/2))]);
end

Alan
Posts: 9536
Joined: December 19th, 2001, 4:01 am
Location: California
Contact:

### Re: Variance Gamma Scaled Self Decomposable process (VGSSD model)

With caveats (looked at the paper only briefly, not a matlab user), it seems from eqn (15) that if you adjust $\theta$ and $\sigma$ for a time step of dt by multiplying them both by $(dt)^{\gamma}$, things might work. I am also assuming your VG simulation code is correct (haven't checked).

I suggest you try it and then, with a large number of simulations, see if all the moments seem to be close to (16)-(19) at various final horizons.

Finally, the values/descriptions of N_sim and T in your posted code looked reversed from what would would expect for that notation.

czar3k
Topic Author
Posts: 11
Joined: October 9th, 2017, 12:34 pm

### Re: Variance Gamma Scaled Self Decomposable process (VGSSD model)

Your remark concernign reversed N_sim & T is right, I was aware of that.
The code itself is correct, I verified it.
Parameters are tuned to daily returns. I set them up as follows:
drift = (theta + mu)*365 = 0.03
sigma*sqrt(365)  = 0.2

Running the code with the following inputs:
VGSimulate(365, 5000, 1, params, S0) produces
366 x 5000 matrix (5000 paths of daily returns)
results in paths that have:
std(log(S(end,:))) = 0.201020089414166
log(mean(S(end,:))) = 0.029067510866851

Swapping input no 1 with niput no 3 that is running the following:
VGSimulate(1, 50000, 365, params, 1)

produces 2 x 50000 matrix of yearly returns

with

std(log(S(end,:))) = 0.200594038070957
log(mean(S(end,:))) = 0.030033611888352

I'm interested in simulation of daily returns, so I have dt = 1.
When I multiply anything  by (dt)^gamma it will remain unchanged.

Alan
Posts: 9536
Joined: December 19th, 2001, 4:01 am
Location: California
Contact:

### Re: Variance Gamma Scaled Self Decomposable process (VGSSD model)

Yes, sorry -- I was way too hasty. Let me take a look at "Self-Decomposability and Option Pricing" and see if I understand it and then get back.

Update. OK, I looked a little more carefully.

These "SSD" processes $X(t)$ are Markov processes with independent increments (pretty sure). But the increments are not (generally) identically distributed and instead depend upon time. This type of process has all the properties of what Sato calls an "additive process".

A simple case would be Brownian motion with time-dependent volatility $\sigma(t)$. If you chose $\sigma^2(t) = (2 \gamma) t^{2 \gamma - 1}$, then you would get $Var\{X(t)\} = t^{2 \gamma}$. Since in that simple case, $X(t) = \int_0^t \sigma(u) dW_u$, that one would be easily simulated. Indeed, one has for all t > s, $X(t) = X(s) + \epsilon_{s,t}$, where $\epsilon_{s,t} \sim N(0, t^{2 \gamma} - s^{2 \gamma})$. I mention this special case because the VGSSD should reduce to this when $\theta = \nu = 0$. Also, it's really the only SSD I think I truly understand at this point.

In general, to do the simulation, you need (in effect) to know the probability transition density $p(t,X(t) | s, X(s))$ and how to (efficiently) make draws from it. For this VGSSD process, I  know neither. If I were doing your project, I would write either Carr, Madan, and/or the authors of the paper you cited with these general questions (not how to fix your matlab code). That later paper seems to be missing a critical appendix cited in footnote 3. Also, it's rather strange that the authors of that one keep referring to "the density function" instead of a whole family of time-inhomogeneous transition densities that seem needed to me.

===============================================================
p.s. If I am reading Sato (Ch. 3) correctly, it seems a brute force simulation would be as follows. Th characteristic function for an additive process increment is

(*)                                        $E[e^{i z (X(t) - X(s))}] = \frac{\phi_{X(t)}(z)}{\phi_{X(s)}(z)}, \quad\quad (t > s)$,

where the numerator and denominator are the VGSSD c.f.'s starting from X(0)=0 that you have in the papers. This follows from Sato Th.16.7 . Now (*) could be Fourier inverted to yield the transition density and drawn from. But this is very general -- holding for all additive processes -- and doesn't use any special VG properties. So, there are likely much more efficient methods. Also, you can check that (*) works for my Brownian example.

czar3k
Topic Author
Posts: 11
Joined: October 9th, 2017, 12:34 pm

### Re: Variance Gamma Scaled Self Decomposable process (VGSSD model)

Thanks for your time on investigation!
I think (just a feeling) that it is impossible that increments (daily returns) are at the same time: independent of previous increments , identically distributed and have a let's call it "kurtosis & skewness" preservation property - central limit theorem forces it. Unless it is a brownian motion as you mentioned  (as it has no excess kurtosis).

So for the moment what I can do is to simulate (independently) a bunch of prices for days 1, 2... 365 so theoretically for any time between 1 and 365 I know the probability of price being above/below certain level and their conditional distribution (conditioned on px level).
Each of those sets of daily, 2-days... 365 days returns will be variance gamma distributed, but  those prices will not be paths...
Unfortunately that's not enough to price barrier options.
I wonder if basing on distribution  of price at time t (say 65 days from now) and distibution of price at t+1d =66 I can workout the distribution of ln(S(66)/S(65)) that is the return itself.
Your remark on probability transition density reminds of local volatility Dupire's model where volatility is non-constant and on step t of simulation depends on S/K and t.

I've already contacted prof. C. O'Sullivan (as I can see on the image in his article there MUST be a way to simulate VGSSD paths) with no reply so far.
I also tried bothering the authors of "Financial modelling - theory implementation and practice"  on Mathworks file exchange . Will try at the source (misters Carr & Madan) , thx for advice.
I wonder why if VGSSD model outperforms VG model as the authors of various articles claim you can not find any code or at least samping procedure description on the internet.
Once again - thanks for the clues - I will keep searching. VGSSD looks very promising in terms of modelling the price evolution especially with MC.

Alan
Posts: 9536
Joined: December 19th, 2001, 4:01 am
Location: California
Contact:

### Re: Variance Gamma Scaled Self Decomposable process (VGSSD model)

czar3k wrote:
Thanks for your time on investigation!
I think (just a feeling) that it is impossible that increments (daily returns) are at the same time: independent of previous increments , identically distributed and have a let's call it "kurtosis & skewness" preservation property - central limit theorem forces it. Unless it is a brownian motion as you mentioned  (as it has no excess kurtosis).

You're welcome.

Yes, i.i.d increments in continuous-time essentially define Levy process. The paper you posted cites another (Konikov & Madan 2002) proving that *all* Levy processes have  c t^(-1/2) and d t^(-1) skewness and excess kurtosis decay. The constants (c,d) are model-dependent and can vanish. Likely true for i.i.d discrete-time as well because (while I haven't actually looked at the cite), the two decay properties likely follow from a short computation with the general form: $E[e^{i z X(t)}]= e^{-t \psi(z)}$.

Of course, the point is that these SSD processes are *not* Levy processes and so need not follow these decay rules.

czar3k
Topic Author
Posts: 11
Joined: October 9th, 2017, 12:34 pm

### Re: Variance Gamma Scaled Self Decomposable process (VGSSD model)

Referring to your brute force simulation idea.
Please let me know I understand correctly.
If the characterictic function of a return of an asset at time t follows VGSSD dist and is expressed as (1-sqrt(-1) *sqrt(t)*nu*theta+0.5*t*nu*sigma^2)^(1/nu) then the characteristic function of dist. of return between t and t+1 (daily return) will be expressed as the ratio of the two c.f., that is:

((1-sqrt(-1) *sqrt(t+1)*nu*theta+0.5*(t+1)*nu*sigma^2)/(1-sqrt(-1) *sqrt(t)*nu*theta+0.5*t*nu*sigma^2))^(1/nu)

I'm  basing my intuition on the example with Gamma dist. from wikipedia on c.f.

If that's correct what transformation do I need to obtain pdf of return from t to t+1?
I'm feeling that we are reaching the boundary of my perception abilities

Alan
Posts: 9536
Joined: December 19th, 2001, 4:01 am
Location: California
Contact:

### Re: Variance Gamma Scaled Self Decomposable process (VGSSD model)

Well, not quite. If X(t) is the mathematical VGSSD process, then

$\phi_{X(t+1) - X(t)}(u)$ is the ratio of the expression on the r.h.s of (15) of O'Sullivan at t+1 and t.

Then, for the risk-neutral stock price process, which has a compensator drift, you need something like

(**) $f_{t,t+1}(u) \equiv \phi_{\log S(t+1) - \log S(t)}(u) = e^{i u r - i u[\omega(t+1) - \omega(t)]} \phi_{X(t+1) - X(t)}(u)$ ,

where $\omega(t) = \log \phi_{X(t)}(u=-i)$.

The pdf of the daily log-return $x = \log S(t+1) - \log S(t)$ is given by the (numerical) Fourier inversion

$p_{t,t+1}(x) = \int e^{-i u x} f_{t,t+1}(u) \frac{du}{2 \pi}$.

Finally, drawing from this pdf can be done by standard "cdf inversion sampling", This requires the cumulative distribution function (cdf), but that too can be obtained by a similar Fourier inversion. I don't have the formulas at hand for you, but you can google for something like "inverting a characteristic function to get a cdf". In this case the c.f. is (**).

Bottom line for the brute force method: set up a function in matlab that produces, via Fourier inversion, the cdf for day t to t+1 returns and sample from it using the cdf inversion sampling method. It won't be fast but it should work and you can then see if the skewness and kurtosis behave as you expect.

p.s. The word "inversion" is used above with two different meanings: "Fourier inversion" and "cdf inversion sampling".  Also, (**) needs to be checked carefully. It looks OK to me now after a half-dozen edits, but no guarantees. It's a daily return version of (3.2) in Carr et al's "Self-Decomposability and Option Pricing"

Alan
Posts: 9536
Joined: December 19th, 2001, 4:01 am
Location: California
Contact:

### Re: Variance Gamma Scaled Self Decomposable process (VGSSD model)

Looking a little more at this interesting process, I think it will take much care to do a correct simulation, no matter what the approach.

If you look at the characteristic function for the VGSSD increments from s to t, you see that it is the ratio of two VGSSD characteristic functions, one using s and one using t. That ratio tends to a constant as $z \rightarrow \infty$. That is a symptom that the distribution of the returns from s to t has a Dirac mass at 0 (assuming $r = \omega = 0$ for simplicity). So, if x is the log-return, the returns density $p_{s,t}(x)$ decomposes into $A \,\delta(x)$ plus a regular part. At a minimum, you need to regularize the characteristic function that is being inverted by subtracting $\phi(\infty)$. This will yield a more suitable function for the integration with some decay as $z \rightarrow \infty$.

You can play around with the special case $\theta = 0, \gamma = 1/2, \nu = 1$, where the inversion is analytically simple and confirm, at least in that case, the Dirac mass.  In that case, $p_{reg,s,t}(x) = a \, e^{- b|x|}$, where $(a,b)$ are easily related to the model parameters and s and t.

czar3k
Topic Author
Posts: 11
Joined: October 9th, 2017, 12:34 pm

### Re: Variance Gamma Scaled Self Decomposable process (VGSSD model)

Thanks for you 2 last post, they inspired me to "play around".
As my math background is not strong enough to deeply feel and understand and rigourously follow your advise (those characteristic functions and numerical integration...) I decided to take different attitude.

As I said I'm in a position to generate yearly returns / prices that are VG distributed I started with doing so. I ended with a vector of 100k final asset prices. Price at time 0 was 1 in each case.
The code is: VGSimulate(1, NPaths, 365, params, 1)
Next step was a little drift from the purely mathematical approach.
I simulated 100k paths of daily returns/prices from the same VG distribution (just scaled drift and sigma accordingly), each path starting not at the price of 1 but at the price taken from previously generated vector of yearly returns/prices.
Next I looked what is the difference between the final price for each of those 100k x 365 days paths and 1 and "corrected" each of 365 prices in those paths to make sure that the final price in all 100k paths is = 1.
Finally I reversed those paths to have S1 = S365, S2 = S364 etc.
I will try with applying correction not to the prices but the returns themselves but that's a different story/
The code was:
multiplier = ([0:365]/365)';

for i = 1:NPaths
[SRev whatever] = VGSimulate(365, 1, 1, params, S(end, i));
delta(i) = 1 - SRev(end);
S2(:,i) = flip(multiplier*delta(i)+SRev);
end

After all I had 100k paths of 365 daily returns which when summed up along the math were VG distributed.
I'm aware that it's not the way I should proceed.
Anyway I checked the properties of the returns for 1st, 2nd,... 365day, they were almost identically distributed (with falling kurtosis as you go from day 1 to day 365, rising std deviation and stable mean). Overall the resuls seem to be quite satisfactory.
As a final check I drew an implied volatility plot which looks good.
As with real IV surfaces vol. falls as time to maturity gets longer and smile effect is stronger for short maturities. For reference I added a flat surface at 0.2 which is the sigma value of VG dist.
Below are some links to the plots so you can have a better picture of what was actually done.
daily paths for T = 365, no correction, yearly returns are normal
daily paths for T = 365 with correction, yearly returns are VG distributed
histogram of D1-D2 returns vs D364-D365 returns, diff. is slight
IV surface

Does it make any sense (I'm sure from the mathematical point of view it's not the right attitude but if it is at least to some extend reasonable with small error margin I'll be satisfied)?

As the daily distribution gradually changes my next idea is to randomly permute daily returns for each path as the order of summing does not affect the sum.
Also I will try to adjust returns instead of prices and see what;s the difference.
I'll post when I'm done it will take a while.

Regards
Thanks for understaning & sorry for my math defficiencies & lack of rigour.
I'm really grateful  for your thourough investigation and help, your post are inspiring even though I miss some points and I'm not able to check your suggestions.

Alan
Posts: 9536
Joined: December 19th, 2001, 4:01 am
Location: California
Contact:

### Re: Variance Gamma Scaled Self Decomposable process (VGSSD model)

Sorry, but I can't follow what you are trying to accomplish here. Whe you say "I checked the properties of the returns for 1st, 2nd,... 365day, they were almost identically distributed (with falling kurtosis as you go from day 1 to day 365, rising std deviation and stable mean). Overall the resuls seem to be quite satisfactory.", this suggests you are just trying to simulate the daily returns of a VG process? But that is easily done by the original code you posted.

On the other hand, if you are trying to approximately simulate the VGSSD process, then the kurtosis shouldn't fall if the approximation is any good.

czar3k
Topic Author
Posts: 11
Joined: October 9th, 2017, 12:34 pm

### Re: Variance Gamma Scaled Self Decomposable process (VGSSD model)

My goal is to simulate VGSSD process at least approximately there is no issue with VG.
Maybe I misunderstand VGSSD process. I think that it is impossible to have everything at once:

1) distributions of daily returns from T(0)->T(1), T(1)->T(2), T(end-1)->T(end) are all skewed with the same skewness
2) distributions of daily returns from T(0)->T(1), T(1)->T(2), T(end-1)->T(end) are all fat tailed with equal kurtosis
3) distributions of daily returns from T(0)->T(1), T(1)->T(2), T(end-1)->T(end) have the same mean and variance
4) distribution of aggregated returns (sum of daily returns from T=0 to N = any value), say early that is from T(0)->T365, have the same skewness & kurtosis as daily returns
5) daily returns are independent
Unless I'm wrong...

Alan
Posts: 9536
Joined: December 19th, 2001, 4:01 am
Location: California
Contact:

### Re: Variance Gamma Scaled Self Decomposable process (VGSSD model)

Yes, I don't think you appreciate the distinction between having independent increments (both VG and VGSSD) and having stationary increments (VG but not VGSSD).

Your properties 1), 2), 3) are false for VGSSD because the daily return distribution changes each day; it is not a stationary increment process. Property 5) is true.

Re property 4), the distribution of the daily increments (daily returns) changes just enough each day to allow what I will call (4)*

(4)* distribution of aggregated returns (sum of daily returns from T=0 to N = any value) have the same skewness & kurtosis as the distribution
of the day 1 returns (the one from T=0 to N=1).

You might want to take a look at Ch, 14 in Cont & Tankov which discusses "additive processes". The VGSSD process is an additive process. For all additive processes, 1), 2), and 3) are generally false;  5) is true.  However, as far as (4)*, while it holds for VGSSD, for other additive models, I think you probably have to do a model-by-model calculation from the characteristic function.

=================================
p.s. Let me repeat my previous p.s. about the characteristic function for these VGSSD increments:

(*)                                 $\phi_{s,t}(z) \equiv E[e^{i z (X(t) - X(s))}] = \frac{\phi_{X(t)}(z)}{\phi_{X(s)}(z)}, \quad\quad (t > s)$,

where the numerator and denominator you have from the linked papers.

So, if you want to see how the various moments of the VGSSD increment distribution are supposed to change from day to day, it is just a tedious exercise in differentiation of the known characteristic function (*). (A good CAS like Mathematica might help, but it is certainly doable by hand). Again from the papers, you already have the moments for the first daily draw: s=0 to t=1. But, I suggest you work out what are the first 4 moments for the daily draw from an arbitrary s to t = s+1. Then, if you are able to develop a good approximate simulation, you can test if your moments change on a daily basis appropriately. Also, I will repeat my caution that the density function for the daily increments has a Dirac mass at zero (with a mass that changes daily), so that is another property I would want to see reproduced in any decent simulator.

czar3k
Topic Author
Posts: 11
Joined: October 9th, 2017, 12:34 pm

### Re: Variance Gamma Scaled Self Decomposable process (VGSSD model)

In that case what would be the expression for pdf of log(S(t+1)/S(t)) for any t in terms of VGSSD characteristic function?
Below if the code for pricing options in VG model, it's ok, I checked it.
I guess the right characteristic function for VGSSD would be:

CF_VG = @(u)(1./(1-1i.*u.*theta.*nu.*sqrt(T) + T.*sigma^2*nu/2.*u.^2)).^(1/nu)

function out= VG_call(type, S0, X, r, T, theta, nu, sigma)

CF_VG = @(u)(1./(1-1i.*u.*theta.*nu + sigma^2*nu/2.*u.^2)).^(T/nu) ;
CF_VG_lgST = @(u) exp(1i.*u.*(log(S0) + r*T + T/nu *log(1 - theta*nu - sigma^2*nu/2))).*CF_VG(u);

integrand1=@(u) real(exp(-1i.* u.* log(X)).* CF_VG_lgST(u-1i)./(1i.*u.* CF_VG_lgST(-1i)));
int11=@(u) real(exp(-1i.* u.* log(X)).* CF_VG_lgST(u-1i)./(1i.*u.* CF_VG_lgST(-1i)));
integrand2=@(u) real(exp(-1i.*u.*log(X)).* CF_VG_lgST(u)./(1i.*u));
int2=@(u) real(exp(-1i.*u.*log(X)).* CF_VG_lgST(u)./(1i.*u));
integral1 = quadgk(integrand1, 0, 10000000 , 'RelTol' ,1e-12, 'AbsTol', 1e-14, 'MaxIntervalCount', 100000);
integral2 = quadgk(integrand2, 0, 10000000, 'RelTol', 1e-12, 'AbsTol', 1e-14, 'MaxIntervalCount', 100000);

Prob1 = 0.5 + (1/pi)* integral1;
Prob2 = 0.5 + (1/pi)* integral2;
call = S0*Prob1 - X*exp(-r*T)*Prob2;
put  = call - S0+exp(-r*T)*X;

if type==1
out = call;
else
out = call - S0+exp(-r*T)*X;
end

Alan
Posts: 9536
Joined: December 19th, 2001, 4:01 am
Location: California
Contact:

### Re: Variance Gamma Scaled Self Decomposable process (VGSSD model)

czar3k wrote:
In that case what would be the expression for pdf of log(S(t+1)/S(t)) for any t in terms of VGSSD characteristic function?

I gave that earlier: (**) on Dec 7 gives the c.f., and right below that is the expression for the pdf (requires a numerical integration)