Serving the Quantitative Finance Community

 
User avatar
Alan
Posts: 2958
Joined: December 19th, 2001, 4:01 am
Location: California
Contact:

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

December 16th, 2017, 5:09 pm

Here is another answer to your original question: Eberlein and Madan give a simulation method in Sec. 5 here (The VGSSD process counts as a Sato process). This is probably your best bet if you can figure out the details (I haven't tried to at this point). But basically, if you know how to simulate a compound Poisson process with drift, you're good to go.
 
czar3k
Topic Author
Posts: 11
Joined: October 9th, 2017, 12:34 pm

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

December 20th, 2017, 11:20 am

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" 

 

 
I think I have some progress. I started with transition from c.f. to PDF of returns at T = 365.
So I have:

CF = 〖(1-iu√t νθ+0.5u^2 tνσ^2)〗^(-1/ν)                                (1)
ω = log⁡(CF(u=-i) ) = log((1-√t νθ-0.5tνσ^2 )^(-1/ν) ) = -1/ν log(1-√t νθ-0.5tνσ^2 )       (2)
PDF= 1/2π ∫〖e^iux∙e^((iur-ω) ) 〗∙CF(u)du                                                                          (3)

Matlab code is as follows:
clear u; clear PDF;


nu = 0.7;
sigma = 0.25/sqrt(365);
T = 365;
X = -0.5:0.005:0.5;
theta = -0.05/sqrt(365);  
r = 0.02/365;


warning off;

for j = 1:numel(X)
    % characteristic function
    CF_T =   @(u) (((1./(1-1i.*u.*theta.*nu.*sqrt(T)+(T).*sigma^2*nu/2.*u.^2)).^(1/nu))) ;
    CF_T2 =  @(u) exp(1i.*u.*(r*T+ 1/nu *log(1 - theta*nu*sqrt(T) - T*sigma^2*nu/2))).*CF_T(u);
    FUN_FINAL = @(u) exp(-1i.*u.*X(j)).* (CF_T2(u));
    % inverse fourier transform
    PDF(j) =  (1/(2*pi)).*real(integral(FUN_FINAL, -Inf, Inf));

end

plot(X, PDF)


It works perfectly.
For the following parameters...

nu = 0.7;

sigma = 0.25/sqrt(365);
T = 365;
X = -0.5:0.005:0.5;
theta = -0.05/sqrt(365);  (needs scaling by sqrt(T) for VGSSD, for VG it scales with T)
r = 0.02/365;

... I received the PDF of yearly returns which you can see here
Here is how the histogram of log(S) generated via  simulation with the code I previously supplied. Match is perfect (I checked for other parameters) so I guess the formulas are free of errors.

Now its time for the second stage which is obtaining PDF of returns from T to T1, at, say T = 1.
The code is as follows, I use "quadgk" instead of "integral" function for speed.
clear u; clear PDF;


nu = 0.7;
sigma = 0.25/sqrt(365);
T = 1;
X = -0.05:0.0005:0.05;
theta = -0.05/sqrt(365);  
r = 0.02/365;


warning off;

for j = 1:numel(X)
    
    CF_T =   @(u) (((1./(1-1i.*u.*theta.*nu.*sqrt(T)+(T).*sigma^2*nu/2.*u.^2)).^(1/nu))) ;
    CF_T_plus_1 =   @(u) (((1./(1-1i.*u.*theta.*nu.*sqrt(T+1)+(T+1).*sigma^2*nu/2.*u.^2)).^(1/nu))) ;
    CF_T2 =  @(u) exp(1i.*u.*(r*(T+1-T) + 1/nu *log(1 - theta*nu*sqrt(T+1) - (T+1)*sigma^2*nu/2) -  1/nu *log(1 - theta*nu*sqrt(T) - (T)*sigma^2*nu/2) )).*(CF_T_plus_1(u)./CF_T(u));
    FUN_FINAL = @(u) exp(-1i.*u.*X(j)).* (CF_T2(u));
    PDF(j) =  (1/(2*pi)).*real(quadgk(FUN_FINAL, -Inf, Inf));
     j
end

plot(X, PDF)
The changes are:

r*(T+1-T) = r    in place of r*T   in  CF_T2

and

+ 1/nu *log(1 - theta*nu*sqrt(T+1) - (T+1)*sigma^2*nu/2) -  1/nu *log(1 - theta*nu*sqrt(T) - (T)*sigma^2*nu/2) 
in place of 
+ 1/nu *log(1 - theta*nu*sqrt(T) - (T)*sigma^2*nu/2)   in  CF_T2

and

*(CF_T_plus_1(u)./CF_T(u))
in place of
*(CF_T(u))   in  CF_T2

Here is what I received, looks strange.
Initially I though that I must have made a mistake in formulas or sth but after playing with integral bounds (I changed them from -Inf to Inf to -500 to 500 which I believe should not change much as PDF values for daily returns beyond this range will be very very close to 0 anyway) I have something like that which looks better and closer to what I expected.
So I guess the problem is with numerical instability of the integral calculated the way as I do it, not with the formula.


I would be grateful for your comments.
 
User avatar
Alan
Posts: 2958
Joined: December 19th, 2001, 4:01 am
Location: California
Contact:

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

December 20th, 2017, 3:26 pm

See my Dec 10 post  
(Suggest you get the easy analytic case given in my post working first; start with pure VGSSD, no stock price drifts).

In the end, though, these one-day increment integrals may be difficult numerically even when regularized. In that case, try the method in my Dec 16 post.
 
czar3k
Topic Author
Posts: 11
Joined: October 9th, 2017, 12:34 pm

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

December 22nd, 2017, 12:02 pm

What exactly is the regularization process about?
I assume it should help to get rid of numerical instability (at least in some cases), so I  guess it's worth trying.
 
User avatar
Alan
Posts: 2958
Joined: December 19th, 2001, 4:01 am
Location: California
Contact:

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

December 22nd, 2017, 2:11 pm

It's about removing the Dirac mass at zero return, which is due to the lack of decay of the Fourier integrand. Your current integrand is not a proper integrand for numerics. Did you read my Dec 10 post?
 
czar3k
Topic Author
Posts: 11
Joined: October 9th, 2017, 12:34 pm

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

December 22nd, 2017, 6:00 pm

Yes, I did, but I have no idea what exactly should be done.
I guess it's about transformation of the integrand function in order to remove the undesired property (Dirac mass) so the integral becomes "numerically processable", Could you provide any example?
 
czar3k
Topic Author
Posts: 11
Joined: October 9th, 2017, 12:34 pm

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

December 22nd, 2017, 6:21 pm

In the meantime I found a more efficient version of matlab code that actually performs inverse Fourier transform on a given function. Article is here
It has some tunable options which I played around with and the results have definetely improved, but the procedure fails for some sets of VGSSD parameters and performs relatively well for the others.
My conclusion is that there is no way to get much better without regularization that you mentioned, because I will be constrained by the memory of my PC or maximum calculation precision which is beyond my control (unless you have a Symbolic Toolbox for Matlab which as far as I know allows for setting any precision). Here are the results
 
User avatar
Alan
Posts: 2958
Joined: December 19th, 2001, 4:01 am
Location: California
Contact:

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

December 22nd, 2017, 7:04 pm

Yes, I did, but I have no idea what exactly should be done.
I guess it's about transformation of the integrand function in order to remove the undesired property (Dirac mass) so the integral becomes "numerically processable", Could you provide any example?
As I suggested a couple times, the case nu=1, theta=0 will illustrate what is going on. The integrand multiplying the exponential is  [$]f(z) \equiv (1 + A z^2)/(1 + B z^2)[$]. To regularize it, you add and subtract the value at z=infinity:

[$] \frac{1 + A z^2}{1 + B z^2} = \frac{A}{B} + \frac{1-A/B}{1 + B z^2}[$]

Now [$]I = \int e^{-i z x} f(z) \frac{dz}{2 \pi}[$] is the sum of two easy integrals, [$]I_1 + I_2[$]. The first one is [$](A/B) \delta(x)[$]. The second one you can do by an easy contour integration (if you have had complex analysis) or, if not, just look it up. It has the exponential form I posted earlier.

Notice that the second integrand is now well-behaved at large z. 

For other (nu,theta), you do the same thing: add and subtract the value at infinity to get a Dirac mass term plus a better behaved integrand for numerical inversion.