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.