Log Normally Scaled Stable Distributions

Needs["StableM7`"];
Needs["StatisticalPlots`"];
data=Uncompress[Import["http://www.mathestate.com/tools/Financial/data/spy.dat","Text"]];
lr=LogReturn[data[[All,2]]];
dtt=Rest[data[[All,1]]];
dtlr=Transpose[{dtt,lr}];
dtlr1=Split[dtlr,Take[#1[[1]],3]===Take[#2[[1]],3]&];
lrp=#[[All,2]]&/@dtlr1;
scl=QkScale/@lrp;

When fitting market data, we noticed a serial dependent structure in the scale factor of intra-day data.  This variation of the scale factor on a daily basis seemed to account for the inability to precisely fit market log return time series to stable distributions.  The plot below shows this variation for SPY on a daily basis starting in July 2007.

DateListPlot[Transpose[{(Last/@dtlr1)[[All,1]],scl}],Joined->True,PlotLabel->"Intra-Day γ",PlotRange->All]

Ignoring the strong serial dependence of the data, we take the sample since the beginning of September 2008, as the parameters have changed since then.

lnpdf=((E^(-((((Log[#1]-#2))^2)/(2 ((#3)^2)))))/((Sqrt[2 π]) #1 #3))&;
scl1=Drop[scl,293];
(*scl1=scl;*)
LNSigma[data_]:=1/1.3489795003921634` Log[Quantile[data,3/4]/Quantile[data,1/4]];
fit=FindMaximum[{Plus@@Log[lnpdf[scl1,μ,σ]],μ∈Reals&&σ∈Reals},{{μ,Median[Log[scl1]]},{σ,LNSigma[scl1]}}];
lndata=Log[scl1];
g1=Histogram[scl1,Automatic,"ProbabilityDensity",PlotRange->All];
g2=Plot[PDF[LogNormalDistribution[fit[[2]][[1]][[2]],fit[[2]][[2]][[2]]],x],{x,0,Max[scl1]},PlotStyle->{Darker[Red],Thickness[0.008]},PlotRange->All];
Show[g1,g2,PlotLabel-> "LogNormal Fit to Intra-Day γ"]
CellPrint[TextCell["MLE fit to lognormal distribution, loglikelihood and parmeters\n"<>ToString[fit],"Text"]];

MLE fit to lognormal distribution, loglikelihood and parmeters
{581.506, {μ -> -7.12784, σ -> 0.505974}}

Below is the quantile-quantile plot of the lognormal fit versus the data.

QQLogNormal[data_,mn_,sd_]:=Module[{nn,qln,ldata},
nn=Length[data];
qln=Quantile[LogNormalDistribution[mn,sd],(Range[nn]-0.5)/nn];
QuantilePlot[qln,scl1,PlotRange->All, AxesLabel-> { "LogNormal","Sample"}]];
Show[QQLogNormal[lndata,fit[[2]][[1]][[2]],fit[[2]][[2]][[2]]],PlotLabel->"Q-Q Plot SPY γ\nLognormal fit"]

This finding suggests that we could model market data as a distribution of the product of a lognormal random variable and a stable random variable.  With that idea we create the characteristic function for just such a mixture of distributions.  This function will have five parameters, the usual α, β, γ, and δ of stable distributions plus σ, representing the standard deviation of the distribution of Log[γ].

Let λ (x, μ, σ) be the lognormal density function.

λ(x,μ,σ)=((e^(-((((log(x)-μ))^2)/(2 (σ^2)))))/((Sqrt[2 π]) x σ))

ϕ(t, α, β) is the standardized stable characteristic function.  (1-parameterization will be used for financial work.)

ϕ(t,α,β)=(e^(-((|t|)^α) (1-i β sgn(t) tan(((π α)/2)))))

A characteristic function for a mixed distribution of characteristic functions is given by (Feller, Vol. 2, p. 504):

(∑_(i=1)^n)(Subscript[p, i]) (Subscript[ω, i])(t)

Where (Subscript[ω, i])(t) is a characteristic function and

(∑_(i=1)^n)(Subscript[p, i])==1

Extending this to a single characteristic function with a varying scale factor, where (Subscript[s, i]) > 0.

(∑_(i=1)^n)(Subscript[p, i]) ω( (Subscript[s, i])t)

If the density function of the scale factor is known and in the case of financial markets may be the lognormal density, we have a market cf (where s > 0; γ is the median of stable γ's).

marketcf(t,α,β,γ,σ)=(∫_0^∞)λ(s,log(γ),σ) ϕ(s t,α,β)ds

Adding a location parameter after the integration, we have a five parameter characteristic function.

marketcf5(t,α,β,γ,δ,σ)=marketcf(t,α,β,γ,σ) (e^(i δ t))

The complete characteristic function is shown below.

MarketCF(t,{α,β,γ,σ,δ})=  ((exp (i t δ))/((Sqrt[2 π]) σ))(∫_0^∞)(1/( s )) exp(-((((log(s)-log(γ)))^2)/(2 (σ^2)))) exp(-((|s t|)^α) (1-(1/2) i s sgn t β ((π α) tan)))ds

This function can be computed by numerical integration and a density from it can be obtained by a double integration of the inverse Fourier transform of the characteristic function.  Such computation is time consuming; a faster approximation of the density can be made using fast Fourier transform of the sampled MarketCF.  It may be possible to set up a series of fast Fourier transforms of stable distributions weighted according to a lognormal distribution, but even this would be slow and cumbersome for maximum likelihood fitting of market data.  It is possible to estimate the parameters for the market distribution by taking advantage of the serial dependence of γ, to approximate γ over many subintervals of a time series.  An example of a fit to market logarithmic returns obtained in this manner is shown below.

LogReturn[data_]:=Differences[Log[N[data]]];
absecf0=Compile[{t,{s,_Real,1}},Abs[Mean[Exp[I t s]]]-1/E];
ABSECF0[t_?NumericQ,s_]:=absecf0[t,s];
QkScale[s_]:= Module[{c0,cc},
c0 =QuartileDeviation[s] ;(*this is best for alpha>1,not a reliable starting point as alpha -> 0*)
(*Absolute value of result is taken in case negative root is found*)
Abs[cc/.FindRoot[ABSECF0[1/cc,s],{cc,c0}]]];
SetOptions[FourierTransform,FourierParameters-> {1,1}];
SetOptions[InverseFourierTransform,FourierParameters-> {1,1}];
SetOptions[Fourier, FourierParameters -> {-1,1}];
SetOptions[InverseFourier, FourierParameters -> {-1,1}];
ReOrder::usage = "ReOrder[list] reorders a frequency domain list for the time domain and vice versa.  If Length[list] is odd the output list will be truncated to the next smallest integer.";
XTable::usage="XTable[SR, N] puts out a list of x values that correspond to the stable x-axis in the time domain for a given length, N, and SR of the DFFT list.";
DFTCF::usage= "DFTCF[cf,SR,N] puts out a discrete FT of a distribution specified by the characteric function cf in the frequency domain compatible with FourierParameters -> {-1, 1}.  SR is the sampling rate; SR is greater than 0.  N is the number of points and must be even.  The output is scaled ready for inversion.";
CF1c::usage= "CF1c[u, α, β, γ, δ] outputs a compiled version of the S(α, β, γ, δ; 1) characteristic function.  All the inputs must be real numbers.";
ReOrder=With[{n=Length[#]},Join[Drop[#,Floor[n/2]],Drop[#,-Floor[n/2]]]]&;
XTable  = Compile[{SR, {N,_Integer}},Module[{n},
n=2 Floor[N/2];
Range[-n/2,n/2-1,1]/(SR)]];
DFTCF [cf_,SR_,n_Integer]:=Module[{rlist,n2},
n2=Floor[n/2];
rlist=2. Pi SR/n Join[Range[0,n2-1],Range[-n2,-1]];
Map[cf[ #]&,rlist]SR/n
	];
CF1c=Compile[{#1,#2,#3,#4,#5},If[#1==0,1,If[#2==1,E^(I*#5*#1-#4*Abs[#1]*((2*I*#3*Log[Abs[#1]]*Sign[#1])/Pi+1)),E^(I*#1*#5-Abs[#1]^#2*#4^#2*(1-I*Sign[#1]*#3*Tan[(Pi*#2)/2]))]]];

data=FinancialData["SPY",All];
pr=data[[All,2]];
lr=LogReturn[pr];
scl2=QkScale/@Partition[lr,50];
cc=Median[scl2];
LNSigma[data_]:=1/1.3489795003921634` Log[Quantile[data,3/4]/Quantile[data,1/4]]
ss=LNSigma[scl2];
parm=SFit[Flatten[Partition[lr,50]/scl2]];
dd=Mean[lr];
cf=DFTCF[Exp[I dd #]Quiet[NIntegrate[((E^(-((((Log[s]-Log[γ]))^2)/(2 (σ^2)))- (s^α)((Abs[#])^α) (1-I β Sign[#] Tan[((π α)/2)])))/((Sqrt[2 π]) s σ))/.{α-> parm[[1]],β->parm[[2]],γ->cc,σ->ss} ,{s,0,∞},Method->{Automatic,"SymbolicProcessing"->0}]]&,1024,2 8192];
pdf=Chop[InverseFourier[cf]];
f=Interpolation[Transpose[{XTable[1024,2 8192],Re[ReOrder[pdf]]}]];
DateListPlot[data,Joined->True,PlotLabel->"SPY Prices"]
g3=Histogram[lr,Automatic,"ProbabilityDensity",PlotRange->All];
g1=Plot[f[x],{x,-0.1,0.1},PlotRange->All,PlotStyle->{Darker[Red],Thickness[0.008]}];
Show[g3,g1,PlotLabel->"Market Function Fit SPY"]
CellPrint[TextCell[ToString[{parm[[1]],parm[[2]],cc,ss,dd}],"Text"]]

{1.86034, -0.0919429, 0.00600552, 0.532775, 0.000232571}

The fit above is the histogram of the log return data with the market function parameters, {α, β, γ, σ, δ}.

The fit below is the stable fit with the stable parameters, {α, β, γ, δ}.  The market function fit refines the differences in the middle of the distribution and has a higher α, that more closely matches the tail behavior of the data.

parm1=SFit[lr];
g2=Plot[SPDF[x,parm1],{x,-0.1,0.1},PerformanceGoal->"Speed",PlotRange->All,PlotStyle->{Darker[Red],Thickness[0.008]}];
Show[g3,g2,PlotLabel-> "Stable Fit"]
CellPrint[TextCell[ToString[parm1],"Text"]]

{1.57356, -0.179032, 0.00604996, -0.0000357222}

We can calculate the log likelihoods for the two fits; the LogNormalStable market function gives a higher log likelihood than the stable fit.

CellPrint[TextCell["Stable log likelihood: "<>ToString[SLogLikelihood[lr,parm1,1]],"Text"]]
CellPrint[TextCell["LogNormalStable log likelihood: "<>ToString[Plus@@Log[f[lr]]],"Text"]]

Stable log likelihood: 12523.7

LogNormalStable log likelihood: 12551.1

txcolor=RGBColor[181/255,18/255,27/255];
t1 = Graphics[Text[Hyperlink[Framed[Style["Financial Home", {txcolor, Bold, Larger}]], "../map/Overview.html"], {-5, 0}]];
t2 = Graphics[Text[Hyperlink[Framed[Style["Next", {txcolor, Bold, Larger}]],  "../s10/LNStableRV.html"], {0, -0.5}]];
t3 = Graphics[Text[Hyperlink[Framed[Style["Download Notebook", {txcolor, Bold, Larger}]], "../s09/LogNormalStable.nb"], {5, 0}]];
t4 = Graphics[Text[Hyperlink[Framed[Style["Previous", {txcolor, Bold, Larger}]], "../wp1/MarketTheory.html"], {0, 0.5}]];

Show[ t1, t2, t3,t4, ImageSize -> 500]
CellPrint[TextCell["\n\n© Copyright " <>ToString[DateList[][[1]]]<>" mathestate    "<>DateString[DateList[],"DateShort"], "Text"]]

© Copyright 2009 mathestate    Fri 16 Jan 2009