Needs["StableM`"];
LogReturn[data_]:=Differences[Log[N[data]]];
AutoCorrelationF[data_,lag_]:=Module[{dataf,dataf1,n},n=Length[data];dataf=Fourier[Join[data-((Plus@@data)/n),Table[0.`,{lag}]]];dataf1=Take[Chop[Fourier[dataf Conjugate[dataf]]],lag+1] N[(Sqrt[n+lag])];(dataf1/(dataf1[[1]]))];
EmpiricCDFList[data_] := Block[{prob,sel,del},
					
					prob= Table[(i-0.5)/Length[data],{i,Length[data]}];
					Transpose[Join[{Sort[data]},{prob}]]];
Clear[DisplayLogStableFit];
scdf[x_?NumericQ,parm_,type_:1,opts___]:=SCDF[x,parm,type,opts];
DisplayLogStableFit[data_,parm_,type_Integer:1]:=Module[{clist,lt,rt,trt,lp1,lp2,lp3,lp4,lp5,nmn,nstd},
clist=EmpiricCDFList[data];
	lt = Select[clist,#[[1]]<0&];
	rt= Select[clist,#[[1]]>0&];
	trt=Transpose[rt];
lp1=ListLogLogPlot[Abs[lt],PlotStyle->Blue,PlotRange->All];
lp2=ListLogLogPlot[Transpose[{trt[[1]],1- trt[[2]]}],PlotStyle->Red,PlotRange->All];
lp3= LogLogPlot[scdf[-x,parm,type],{x,0.0001,-4 Min[data]},PlotStyle->Blue,PlotRange->All];
lp4= LogLogPlot[1-scdf[x,parm,type],{x,0.0001,-4 Min[data]},PlotStyle->Red,PlotRange->All];
nmn=Mean[data];
nstd=StandardDeviation[data];
lp5= LogLogPlot[1-scdf[x,{2,0,nstd /Sqrt[2],nmn},1],{x,0.0001,4 nstd},PlotStyle->Green,PlotRange->All];
Show[lp1,lp2,lp3,lp4,lp5,PlotLabel->"Log Log Left Tail Blue, Right Red, (Normal Green)",AspectRatio->Automatic,PlotRange->All]];
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}]]];
Clear[EstimateScale];
EstimateScale[s_List,α_,β_]:=MeanDeviation[s]/(Gamma[1-1/α] Cos[ArcTan[β Tan[α Pi/2]]/α]/(Pi/2  Cos[ArcTan[β Tan[α Pi/2]]]^(1/α)));
ScaleS1[n_,s1params_List] := {s1params[[1]],s1params[[2]],n^(1/s1params[[1]])s1params[[3]],n s1params[[4]]};

Volatility Model

Our theory of markets is derived from the continuous double auction model (CDA), with the requirement that the buy and sell limit order books organize the prices in a manner that each has a heavy power tail.  With the CDA model, price movements occur when there is a mismatch between buy and sell market orders; the excess orders are then filled from the appropriate limit order book.  If the limit order books have price differences with spacing that will yield heavy tailed returns when executed, and if the events in the market place are independent, then by the generalized central limit theorem (GCLT) sums of the returns will converge to a stable distribution.  We have further shown that a pure power tail random variable can arise from taking the exponent of an exponential random variable.  For all this to be useful, we need a few more assumptions.  First it would be nice if the exponent of the power tail distribution in the order books remained constant or nearly constant for a given security over time.  It would be speed convergence if the two tail exponents were equal, but this is not a necessary feature for stable distributions to arise from the model.  There must also be enough transactions for the sums to begin to converge to a stable regime; how many transactions are needed is uncertain for it seems that there is convergence to a distribution that we cannot distinguish from stable with relatively few sums.  To get a stable distribution with the same tail exponent as the distributions from which the sums arise, does take a large number of sums.

For the model below we assume that the tail exponent is constant and we further assume that β is determined by the relative ratios of how many unmatched market orders are filled in each of the limit order books.  We assume this to be constant, but it need not be.  With all the assumptions above sums of returns will by the GCLT converge to stable distributions.

On the Market Data Fitting page, we proposed that actual stock market returns could be modeled by a non-stationary stable distribution model where the scale factor was continually varying, with patterns that showed serial dependence.  We now study that idea with pseudo-random stable data.  We will take a large random sample generated with known parameters and modulated with known functions; then see how well we can recover the original parameters.

Using the SPY intraday minute by minute data we saw that there was a pattern to the intraday trading.  The very first return of the day is actually not a minute return, but the return from the previous close, so it has a much larger scale factor than the minute by minute returns during the day.  The rest of these returns tend to a cycle of higher volatility at the beginning and end of the day and a low point mid day.  For our example we use a parabolic function.

idcycle=ReplacePart[(1.6061696252465485`-0.006185404339250495` #1+0.00001577909270216963` ((#1)^2))&/@Range[391],1-> 10];
intradaycycle=idcycle/Mean[idcycle];
ListPlot[intradaycycle,PlotRange->All]

We have purposely adjusted this intraday scale factor so it will have a mean of one so we can show graphs of the recovered function without rescaling.  The first return has a very large scale factor, but all the other minute returns throughout the day have a smoothly changing scale factor.

Mean[intradaycycle];

Below we create another function to model the variation of day to day volatility.  We use summed trigonometric functions for our known function.  This is smoothly varying as we might expect for day to day volatility in the market place.  We choose a daily scale factor from this function.

volfn=1.`+0.125` Cos[((2 π #1)/63)]+0.0625` Cos[((2 π #1)/21)]+0.5` Sin[((π #1)/126)] &;
Plot[volfn[x],{x,1,252}]

vol=volfn/@Range[252];
ListPlot[vol]

We have again set the mean of this function to one, so we will be multiplying by functions with an average scale of one.  This is not needed for the model, but is done purposely because our scaling function adjusts to a stable scale factor of one and we want to test its performance.

Mean[vol];

Next we create a stable random sample of size 98,532 corresponding to 252 trading days each containing 390 minute intervals plus 1 interday interval.  The starting parameters are {α, β, γ, δ} = {1.8, -0.1, 1, 0}.

data=SRandom[252 391,{1.8,-0.1}];

data1=Partition[data,391];

To give the final distribution parameters similar to the real life SPY data we multiply by the factor below.

cq=0.006 391^(-1/1.8)

The plot below shows the data rescaled with the varying scale model.  The intraday cycle is imperceptible, but the higher volatility of the daily volatility function can be seen at the beginning of the series, volatility then falls and rises at the end of the series--very much like real market data.

data2=cq Flatten[intradaycycle #&/@(vol data1)];
ListPlot[data2,Joined->True,PlotRange->All,PlotLabel->"Simulated Market Returns"]

The autocorrelation plot of the absolute value of each data point shows the serial dependent structure that we have created along with the intraday cycle and the spike between each day's cycle.

ListPlot[AutoCorrelationF[Abs[data2],391 50],Joined->True,PlotRange->All,PlotLabel-> "Autocorrelation Function Abs[LogReturn] - 50 Day Lag"]

If we fit this data sample as if it were drawn from a stationary stable distribution, we find exactly the same phenomenon that we observed with actual market data.  Alpha is artificially low and the plot of the tail fit shows that the data tails are lighter than would be predicted by a stable model.

p1=Quiet[SFit[data2]];
CellPrint[TextCell["Stable parameters\n"<>ToString[FortranForm[#]&/@(Round[p1 10000.]/10000)],"Text"]];
Quiet[DisplayLogStableFit[data2, p1]]

Stable parameters
{1.5879, -0.0484, 0.0002, 0}

Now we use our scaling method to recover the volatility functions.  First we go after the day to day data.  The plot shows the scale factor data points we get for each point along with our known volatility function  plotted in red.

scl=QkScale/@Partition[data2,391];
g1=ListPlot[scl];
g2=Plot[volfn[t] cq,{t,1,252},PlotStyle->Red];
Show[g1,g2]

But since we know the function with which we started, we can test our accuracy.  We look at the differences between our data points and find that they are approximately normally distributed and have a standard deviation that is about 5% of our mean scaling function

diff=vol-scl/cq;
CellPrint[TextCell["Standard Deviation Errors: "<>ToString[StandardDeviation[diff]],"Text"]];
CellPrint[TextCell["Mean Scaled Errors: "<>ToString[Mean[scl/cq]],"Text"]];
CellPrint[TextCell["Mean Volatility Function: "<>ToString[Mean[vol]],"Text"]];

Standard Deviation Errors: 0.04846370661524439

Mean Scaled Errors: 0.9714122058379646

Mean Volatility Function: 1.

Now we look at the autocorrelation function with the day to day volatility removed and find that we have left the intraday cycle.

sdata=Flatten[Partition[data2,391]/scl];
ListPlot[AutoCorrelationF[Abs[sdata],391 15],Joined->True,PlotRange->All,PlotLabel->"Autocorrelation Function Abs[LogReturn] - 15 Day Lag"]

Next we recover the intraday scale function.  The known function is in red and the recovered points in blue.

scl2=Quiet[QkScale/@Transpose[Partition[sdata,391]]];
g1=ListPlot[intradaycycle,Joined-> True,PlotStyle->{Red,Thickness[0.005]},PlotRange->All];
g2=ListPlot[scl2,PlotStyle->Blue,PlotRange->All];
Show[g2,g1]

sdata1=Flatten[Transpose[Transpose[Partition[sdata,391]]/scl2]];
ListPlot[sdata1,Joined->True,PlotRange->All,PlotLabel-> "LogReturns with Serial Dependence Structure Removed"]

The fully rescaled data is shown above and the autocorrelation function of the absolute value of the simulated returns is below, no trace of the serial dependence functions remain.

ListPlot[AutoCorrelationF[Abs[sdata1],391 15],Joined->True,PlotRange->All]

We fit the data and recover the original α and β with surprising accuracy.  The plot of the fit shows true stable tails.

p2=SFit[sdata1];
CellPrint[TextCell["Stable parameters\n"<>ToString[FortranForm[#]&/@(Round[p2 10000.]/10000)],"Text"]];
Quiet[DisplayLogStableFit[sdata1,p2]]

Stable parameters
{1.8094, -0.1237, 1.003, -0.0034}

Here is the fit to the original data before the volatility model was superimposed on it.  The result is virtually the same.

p3=SFit[data];
CellPrint[TextCell["Stable parameters\n"<>ToString[FortranForm[#]&/@(Round[p3 10000.]/10000)],"Text"]];
Quiet[DisplayLogStableFit[data,p3]]

Stable parameters
{1.8124, -0.1327, 1.0024, -0.0049}

We took a stable random sample containing nearly a hundred thousand points and multiplied every point in the sample by a serially dependent structure.  With knowledge of stable scaling functions we were able to accurately recover the original parameters.  In the model we have ignored the stable δ parameter; there are several reasons for this.  On a minute by minute basis the value will be so small that we cannot likely detect it accurately.  While the market system might build in an expectation of inflation and inflation adjusted return, we suspect that it would always be changing.  Stock prices seem to be adjusted for these factors when the companies report their quarterly earnings.  Thus long range expectations need not be included in the market model.

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}]],  "../r08/Crash.html"], {0, -0.5}]];
t3 = Graphics[Text[Hyperlink[Framed[Style["Download Notebook", {txcolor, Bold, Larger}]], "../r07c/VolatilityModel.nb"], {5, 0}]];
t4 = Graphics[Text[Hyperlink[Framed[Style["Previous", {txcolor, Bold, Larger}]],  "../r07/Volatility.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 2008 mathestate    Fri 22 Feb 2008