Testing The Model

Our central idea is that market price formation is driven by sums of logarithmic returns of many transactions.  We believe that the price formation process of the continuous double auction has heavy tails.  Over times as short as a minute, the returns begin to converge to a stable distribution by the generalized central limit theorem.  When we try to fit long stretches of returns over many days, we quickly find that the data do not conform to a stationary stable distribution.  The tails of the return data are too light.  We either have the wrong distribution or bad methods.  By experimentation we find that the method of maximum likelihood fitting of the distribution shape parameter, α, can be significantly affected by contamination with an excess of returns in the middle of the distribution.  This distorts the calculation of α, even though the tails of the contaminated data are not changed.  When we look at events that are in the process of converging to a stable regime, by a two tailed continuous double auction model, we find that there are initially an excess of zero returns in the model, upon summation these disappear, but to the extent that a sample has not fully converged to a stable regime, an excess of small returns artificially lowers α.

The experiment below shows a stable random sample generated with parameters: {1.8, 0.2, 1,0 }. It is fit by the maximum likelihood method below.

Needs["StableM7`"];
(*Needs["Histograms`"]; Version6.0*)
s=SRandom[1000,{1.8,0.2,1,0}];
pp1=SFit[s];
CellPrint[TextCell["Stable Parameters {α, β, γ, δ} = "<>ToString[pp1],"Text"]];

Stable Parameters {α, β, γ, δ} = {1.80643, 0.535578, 1.03628, 0.0723559}

The sample is contaminated with a hundred uniform random variables range {-0.2, 0.2} and refit.  The graph below shows the plot of the density of the two fits.  Calculated α is lower in the sample with extra mass added to the center of the distribution, even though the mass on the tails has not changed.

rr=RandomReal[{-0.2,0.2},100];
pp2 =SFit[Join[s,rr]];
CellPrint[TextCell["Stable Parameters {α, β, γ, δ} = "<>ToString[pp2],"Text"]];
Plot[{SPDF[x,pp1],SPDF[x,pp2]},{x,-10,10},PlotRange->All,PlotStyle->Thickness[0.003], PlotLabel-> "Blue - Pure Stable Sample\nRed - Contaminated Sample"]

Stable Parameters {α, β, γ, δ} = {1.77813, 0.438959, 0.965275, 0.0620158}

In high resolution minute by minute data, this phenomenon may be contributing because the sums of log returns have not fully converged to a stable regime.  We believe, there is more to the problem, namely that the scale factor, γ, of the distribution is changing on a minute by minute basis.  The ultimate collection of returns we see is derived from a mixture of stable distributions with varying γ.  The variation of the scale parameter, γ, is not random, but shows serial dependence with an interesting pattern.  The blue curve is the autocorrelation function of the raw absolute value of the log returns.  Our data set consists of 391 returns a day minute by minute.  It does not include the pre-market or after-market trading.  The log return from the previous close to the open is included as the first return of the day and all the days are concatenated.  The spikes are caused by the inter-day returns, but there is an intra-day pattern that is caused by higher volatility at the begin and end of the day.  From day to day there is serial dependence that slowly decays.  Fifteen days of data are shown in the graph.  The red graph divides each day's returns by the scale factor for that day.  There is no attempt to correct the cyclic intra-day variation.  With this maneuver, most of the day to day serial dependence disappears, but the intra-day cycle and inter-day spikes remain.

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]]))];
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}]]];
symbol="SPY";
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;
lr1=Flatten[lrp/scl];
ListPlot[{AutoCorrelationF[Abs[lr],15 391
],AutoCorrelationF[Abs[lr1],15 391
]},PlotRange->All,Joined->True,PlotLabel->"Blue - Minute Data Autocorrelation Abs[LogReturn]\nRed -  Scale Adjusted",Frame->True,Axes->False]

The next graph shows the plot of the SPY ETF closing price in blue and the daily stable γ, calculated since July 2007.  The spike in volatility at the end of the series is dramatic and historical.  There has not been such a spike in the history of the SPY ETF and a look at other historical data such as the Dow Jones Industrial Average, suggests this level of volatility has not occurred since the 1930s.  As this page is updated, it will be interesting to see how this plays out.  Right now it presents an opportunity to test our model.  Note that spikes of volatility (γ) are associated with price declines.  Periods of falling volatility are associated with rising prices.

dtpr=Last/@Split[Rest[data][[All,{1,2}]],Take[#1[[1]],3]===Take[#2[[1]],3]&];
dt=Last[#][[1]]&/@dtlr1;
prticks=Table[{i,2.5i+80},{i,0,88,4}];
gticks=Table[{i,i/10000.},{i,0,30,5}];
DateListPlot[{{#1,(#2-80)/2.5}&@@@dtpr,Transpose[{dt,10000scl}]},Joined->True,PlotLabel->symbol<>" Close/Gamma",FrameTicks-> {Automatic,prticks,Automatic,gticks},PlotRange->All]

Our hope is that our relatively simple model of varying scale factor will have a relatively constant shape parameter, α.  Such a finding would make analysis of the problem somewhat mathematically tractable.  We expect the β parameter to vary on a minute to minute basis dependent upon the ratio of orders filled in the buy and sell order books.  On a minute to minute basis, δ should be essentially 0.

Needs["Stable`"];
(*Stable.m may be obtained from www.robustanalysis.com*)
parm=Quiet[StableMLFit/@lrp];
gq0=DateListPlot[Transpose[{dt,parm[[All,1]]}],Joined->True,PlotRange->All,PlotLabel->symbol<>" Alpha"];
{ma,mb,mc,md}= Mean[parm];
ci=StableMLEConfidenceInterval[{ma,mb,mc,md},391];
gq0a=Plot[{ma+ci[[1]],ma-ci[[1]]},{x,AbsoluteTime[First[dt]],AbsoluteTime[Last[dt]]},PlotStyle-> {{Thickness[0.008],Darker[Red]}}];
CellPrint[TextCell["Mean stable parameters "<>ToString[AccountingForm[{ma,mb,mc,md},NumberSigns->{"-", " "}]],"Text"]];
Show[gq0,gq0a]

Mean stable parameters { 1.82031, -0.0375413,  0.000477175, -0.00000580621}

The graph above is the calculated value of α for each day of the sample.  The heavy red lines are the 95% confidence intervals based on the mean α and a sample size for each day of 391.  Some days are actually shorter because of early market closing before holidays.  There are more excursions outside the confidence intervals than we would predict, but it is encouraging that most of the time α is close to the mean, especially at the end of the series where we have begun to experience extreme volatility.  The scale factor is changing, but the basic fractal dimension, α, of the market behavior seems to be rather constant.  Some of the variation and low values for α calculated may have been from excess zero returns in the sample.  Note that there are fewer zero returns at the end of the series when volatility and market volume were high, suggesting more complete convergence to a stable regime at the one minute time interval.

DateListPlot[Transpose[{dt,Length[Select[#,#==0&]]&/@lrp}],PlotLabel-> "Number of Zero Returns Each Day"]

We can use extreme value distributions and theory to make an estimate of α from the tail data, thereby overcoming some of the problem of too much mass in the middle of the mixture sample of our data.  The idea here is that tail events occur with some regularity.  So we have selected the maximum values from each day's morning and afternoon events and fit this distribution to a generalized extreme value distribution (gev).  The ξ parameter of the gev is equivalent to 1/α of a stable distribution from which it is derived; thus we have another method of calculating α.

Before September 2008, this calculation worked well for the whole database, with the surge of volatility since then the magnititude extreme events in the last part of the data set overwhelms the calculculation.  We believe the distribution of the intra-day gamma has significantly changed.  This calculation is now being performed on the last part of the data that begins with September 2008.

Clear[loggev];
loggev[x_,ξ_,c_,d_]:=-(((1+(((-d+x) ξ)/c)))^(-1/ξ))-Log[c]-(((1+ξ) Log[1+(((-d+x) ξ)/c)])/ξ)/;ξ≠0
loggev[x_,ξ_,c_,d_]:=-(E^((d-x)/c))+((d-x)/c)-Log[c]/;ξ== 0
SetAttributes[loggev,Listable]
gev[x_,ξ_,σ_,μ_]:=(Piecewise[{{(((E^(-(((1+(((x-μ) ξ)/σ)))^(-1/ξ)))) (((1+(((x-μ) ξ)/σ)))^(-1-(1/ξ))))/σ), ξ≠0&&x ξ-μ ξ+σ>0}, {((E^(-(E^((-x+μ)/σ))+((-x+μ)/σ)))/σ), ξ==0&&σ>0}}]);
ss2=Flatten[{Max[Take[#,Floor[Length[#]/2]]],Max[Take[#,-Length[#]+Floor[Length[#]/2]]]}&/@Drop[lrp,293]];
ulim=Max[ss2];
llim=Min[ss2];
w=Quiet[FindMaximum[{Plus@@loggev[ss2,a,c,d],c+a llim>a d&&c+a ulim>a d},{{a,1/2},{c,0.003},{d,0.001}}]];
{aa,cc,dd}={a,c,d}/.w[[2]];CellPrint[TextCell["Sample size of right tail "<>ToString[Length[ss2]]<>".","Text"]];
h1=Histogram[ss2,Automatic,"ProbabilityDensity"];
g1=Plot[gev[x,aa,cc,dd],{x,llim,ulim},PlotStyle->{Thickness[0.008],Darker[Red]},PlotRange->All];
Show[h1,g1, PlotLabel->"GEV Fit to Right Tail Data"]
CellPrint[TextCell["Calculated stable α "<>ToString[1/aa]<>".","Text"]];

Sample size of right tail 210.

Calculated stable α 1.9203.

ss3=Flatten[{-Min[Take[#,Floor[Length[#]/2]]],-Min[Take[#,-Length[#]+Floor[Length[#]/2]]]}&/@Drop[lrp,293]];

ulim=Max[ss3];
llim=Min[ss3];
w=Quiet[FindMaximum[{Plus@@loggev[ss3,a,c,d],c+a llim>a d&&c+a ulim>a d},{{a,1/2},{c,0.003},{d,0.001}}]];
{aa,cc,dd}={a,c,d}/.w[[2]];
CellPrint[TextCell["Sample size of left tail "<>ToString[Length[ss2]]<>".","Text"]];
h1=Histogram[ss3,Automatic,"ProbabilityDensity"];
g1=Plot[gev[x,aa,cc,dd],{x,llim,ulim},PlotStyle->{Thickness[0.008],Darker[Red]},PlotRange->All];
Show[h1,g1, PlotLabel->"GEV Fit to Left Tail Data"]
CellPrint[TextCell["Calculated stable α "<>ToString[1/aa]<>".","Text"]];

Sample size of left tail 210.

Calculated stable α 1.57121.

At this point with a smaller data set the tails are different.  We expect that as we accumulate more data at the current volatility levels they will become closer.  Before the market volatility in September 2008, these fits seemed more consistent with a stationary α at approximately 1.8.  

When we fit the whole raw data set, we get the following result, with the finding of an artificially very low α.  The scaled data find an α more consistent with the daily fits, but now the γ and δ are scaled to γ ≅ 1.

pp0=StableMLFit[lr];
CellPrint[TextCell["Stable fit to raw data {α, β, γ, δ} = "<>ToString[AccountingForm[pp0,NumberSigns->{"-", " "}]],"Text"]];
pp1=StableMLFit[Flatten[lrp/scl]];
CellPrint[TextCell["Stable fit to scaled data {α, β, γ, δ} = "<>ToString[AccountingForm[pp1,NumberSigns->{"-", " "}]],"Text"]];

Stable fit to raw data {α, β, γ, δ} = { 1.38814,  0.00000000280358,  0.000367778, -0.0000000972526}

Stable fit to scaled data {α, β, γ, δ} = { 1.79855,  0.00000000988973,  1.01309,  0.000111789}

We now analyze the properties of the scale factor, γ, which is a measure of market volatility.  The graph below shows the plot of γ, calculated by two methods one is a quick method based on the stable characteristic function the other is by the maximum likelihood method--the results are virtually identical.

gq1=ListPlot[{scl,Transpose[parm][[3]]},Joined->True,PlotRange->All,PlotLabel->symbol<>" Gamma\nBlue - Quick Routine\nRed - ML Fit"]

The pattern of volatility is related to market volume, but they are not directly proportional .  The following plot shows the daily volume over the same time interval as the sample.

vol=FinancialData[symbol,"Volume",First[data][[1]]];
DateListPlot[vol,Joined->True,PlotRange->All,PlotStyle->Black,PlotLabel-> symbol<>" Volume"]

We show below the histogram of the daily γ data beginning September 2008.  These data are clearly not random.  The idea of showing their histogram and a fit to a lognormal distribution is mainly to show that the data are somewhat constrained.  Until the most recent week, the sample was well fit by the lognormal distribution, but the extreme results of the recent week are too heavy to be consistent with the rest of the sample and a lognormal fit.  (The data before September 2008, were also fit with a log normal distribution, but the parameters have changed, so we are tracking the more recent data set.)

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"];
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
{678.354, {μ -> -7.15085, σ -> 0.482587}}

The Q-Q plot is shown.  The lognormal fit was quite good until the current weeks, when the divergence of the tail began with the current high volatility.  We are likely in a new era of volatility which will persist for some time.  The intra-day scale factor may not be sufficiently constrained to fit to a model such as a lognormally scaled stable distribution.

Needs["StatisticalPlots`"];
QQLogNormal[data_,mn_,sd_]:=Module[{nn,qln,ldata},
nn=Length[data];
qln=Quantile[LogNormalDistribution[mn,sd],(Range[nn]-0.5)/nn];
QuantilePlot[qln,data,PlotRange->All, AxesLabel-> { "LogNormal","Sample"}]];
fit1=FindMaximum[{Plus@@Log[lnpdf[scl1,μ,σ]],μ∈Reals&&σ∈Reals},{{μ,Median[Log[scl1]]},{σ,LNSigma[scl1]}}];
Show[QQLogNormal[scl1,fit[[2]][[1]][[2]],fit[[2]][[2]][[2]]],PlotLabel->"Q-Q Plot SPY γ\nLognormal fit"]

We are quite pleased with the performance of the model over this very extreme market period, but the constraint of a lognormal distribution for the scale factor seems to have different parameters since September 2008.  We will update this page periodically.  This model simply is a non-stationary stable distribution with a varying scale factor or amplitude.  Since stable distributions can have heavy tails they are prone to extreme jumps.  The varying scale factor is multiplicative so the jumps can be magnified greatly. Over the last year we have seen this scaling increase by a factor of 13.  The scaling always in the past has shown serially dependence, so the current volatility is not likely to disappear soon.  By this model extreme events are common when volatility is high and crashes should be expected.

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["Market Theory", {txcolor, Bold, Larger}]],  "../wp1/MarketTheory.html"], {0, 0}]];
t3 = Graphics[Text[Hyperlink[Framed[Style["Download Notebook", {txcolor, Bold, Larger}]], "../s13/ModelTest.nb"], {5, 0}]];
(*t4 = Graphics[Text[Hyperlink[Framed[Style["Previous", {txcolor, Bold, Larger}]], "../s10/LNStableRV.html"], {0, 0.5}]];*)

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

© Copyright 2009 mathestate    Sun 1 Feb 2009