Fitting Market Data

Needs["Stable`"];
LogReturn[data_]:=Differences[Log[N[data]]];
ScaleS1[n_,s1params_List] := {s1params[[1]],s1params[[2]],n^(1/s1params[[1]])s1params[[3]],n s1params[[4]]};
EmpiricCDFList[data_] := Block[{prob,sel,del},
					
					prob= Table[(i-0.5)/Length[data],{i,Length[data]}];
					Transpose[Join[{Sort[data]},{prob}]]];

ft={0.5,0.1,0.05,0.01,0.005,0.001,0.0001};
lfticks=Transpose[{ft,ft}];
rfticks=Transpose[{ft,1-ft}];
Clear[DisplayLogStableFit];
DisplayLogStableFit[data_,parm_,type_Integer:1]:=Module[{clist,lt,rt,trt,lp1,lp2,lp3,lp4,lp5,nmn,nstd,xlim,xlim0},
clist=EmpiricCDFList[data];
	lt = Select[clist,#[[1]]<0&];
	rt= Select[clist,#[[1]]>0&];
	trt=Transpose[rt];
xlim0=Min[rt,Abs[lt]];
xlim=Max[{-Min[data],Max[data]}];
lp1=ListLogLogPlot[Abs[lt],PlotStyle->Blue,PlotRange->All,Frame->True,FrameTicks->{True,lfticks,None,rfticks},Axes->False];
lp2=ListLogLogPlot[Transpose[{trt[[1]],1- trt[[2]]}],PlotStyle->Red,PlotRange->All,Frame->True];
lp3= LogLogPlot[StableCDF[-x,parm,type],{x,xlim0,xlim},PlotStyle->Blue,PlotRange->All];
lp4= LogLogPlot[1-StableCDF[x,parm,type],{x,xlim0,xlim},PlotStyle->Red,PlotRange->All];
nmn=Mean[data];
nstd=StandardDeviation[data];
lp5= LogLogPlot[1-StableCDF[x,{2,0,nstd /Sqrt[2],nmn},1],{x,xlim0,4.5 nstd},PlotStyle->Green,PlotRange->All];
Show[lp1,lp2,lp3,lp4,lp5,PlotLabel->"Log Log Stable Distribution\nLeft Tail Blue, Right Red, (Normal Green)",PlotRange->All,FrameLabel->{"Abs[x]","Left Tail",None,"Right Tail"}]];

A reasonably accurate fit of stable parameters to a data set can be obtained by the maximum likelihood method.  Below are plots of the SPY prices since 1993, the calculated density function with histogram, and the empirical density and distribution.  The calculated fits are in blue.

data=FinancialData["SPY",All];
DateListPlot[data,PlotLabel->"SPY Prices",Frame->True,Joined->True]
lr=Differences[Log[data[[All,2]]]];
g3=Histogram[lr,Automatic,"ProbabilityDensity",PlotLabel->"SPY Daily Log Returns",ChartStyle->RGBColor[0.5`,0.5`,0.5`],PlotRange->All,Frame->True];
pr=(Range[Length[lr]]-0.5)/Length[lr];
{First[pr],Last[pr]};
g4=ListPlot[Transpose[{Sort[lr],pr}],PlotLabel->"Empirical Probability Distribution",Frame->True,PlotStyle-> Gray];
parm = StableMLFit[lr,1];
(*TableForm[{parm},TableHeadings-> {None,{"α","β", "γ","δ"}}]*)
g7=Plot[StablePDF[x,parm,1],{x,Min[lr],Max[lr]},PerformanceGoal->"Speed",PlotStyle->{Darker[Red],Thickness[0.005]},PlotRange->All,Frame->True];
g9=Show[g3,g7,PlotLabel->"Stable Density Fit",Frame->True];
g8=Plot[StableCDF[x,parm,1],{x,Min[lr],Max[lr]},PerformanceGoal->"Speed",PlotStyle->{Darker[Red],Thickness[0.005]},PlotRange->All,Frame->True];
g10=Show[g4,g8,PlotLabel->"Stable Distribution Fit"];
Show[GraphicsRow[{g9,g10}]]

The fit of financial data to stable distributions has some problems.  Above we show some curves of the fit of SPY data to stable density and distribution curves.  Those pictures, however, don't well show the detail of the tails.  Below are the stable parameter fit using maximum likelihood method.  The plot beneath that shows a log-log fit to the cumulative distribution function; to align the tails, the absolute value of the left tail is shown in blue; the right tail is shown as 1- probability in red.  In this format stable distributions with α < 2, show linear parallel tails with the slope of minus α.  If β is zero the tails are superimposed (when β = ±1, the lighter tail is not linear).  The tail of a Normal distribution calculated from the sample is shown in green and is not linear.  The dots show the data points of the respective tails.

TableForm[{parm},TableHeadings-> {None,{"α","β","γ","δ"}},TableSpacing-> {1,5}]
Quiet[DisplayLogStableFit[lr,parm,1]]

The actual data tails are too light with slope steeper than that of the fit.  When the actual slopes of the tails of the data are calculated from the data tails (upper and lower 2%), the tail exponents are far too high (graph below).  To fit stable laws tail exponents must be less than 2.  By the generalized central limit theorem sums of random variables from distributions with linear tail exponents that are less than two will converge to a stable distribution.  Those with higher tail exponents will converge to a Normal distribution on  summation.  Data with tail exponent higher than two cannot be from a stable distribution.

cdf=EmpiricCDFList[lr];
qp1=ListPlot[{Log[{Abs[#1],#2}]&@@@Take[cdf,Floor[Length[cdf]/50]],Log[{Abs[#1],1-#2}]&@@@Take[cdf,-Floor[Length[cdf]/50]]},PlotStyle-> {Blue,Red}];
yleft=Fit[Log[{Abs[#1],#2}]&@@@Take[cdf,Floor[Length[cdf]/50]],{t,1},t];
yright=Fit[Log[{Abs[#1],1-#2}]&@@@Take[cdf,-Floor[Length[cdf]/50]],{t,1},t];
qp2=Plot[{yleft,yright},{t,-4,-2},PlotStyle-> {Blue,Red}];
Show[qp1,qp2,PlotLabel-> "Right Red, α = "<>ToString[-yright[[2]][[1]]]<>"\nLeft Blue, α = "<>ToString[-yleft[[2]][[1]]],AxesLabel-> {"Log[LogReturn]","Log[Probability]"}]

Conversely, if we take a simulation derived from the fit of financial data to a stable distribution with parameters such as above and reconstruct a price series, we find extreme changes more frequently than are ever found in a financial time series.  Thus the stable fit over estimates the frequency of extreme events, for financial data.  There is a dynamic demonstration of the graph below on the Wolfram Demonstrations Project.

ListPlot[100 Exp[FoldList[Plus,0,StableRandom[1000,parm,1]]],Joined ->True,PlotLabel-> "Daily Price Simulation"]

Another problem in using a statistical distribution for financial data is that financial data are not completely random.  Below an autocorrelation plot of the SPY price series of log returns looks random, i.e. no observable autocorrelation.

AutoCorrelationF[data_,lag_]:=Module[{dataf,dataf1},dataf=Fourier[Join[data-Mean[data] ,Table[0.`,{lag}]]];
dataf1=Take[Chop[Fourier[Abs[dataf]^2]],lag+1] ;
(dataf1/(dataf1[[1]]))];ListPlot[AutoCorrelationF[lr,126],PlotRange->All,Joined->True,PlotLabel->"One Day Autocorrelation LogReturn",Frame->True,Axes->False]

However, when you look at the same autocorrelation function with the absolute value of the log returns there is evidence of serial dependence in the data.

ListPlot[AutoCorrelationF[Abs[lr],126],PlotRange->All,Joined->True,PlotLabel->"One Day Autocorrelation Abs[LogReturn]",Frame->True,Axes->False]

If you scramble the order of the same data above the serial dependent structure disappears.

ListPlot[AutoCorrelationF[Abs[RandomSample[lr,Length[lr]]],126],PlotRange->All,Joined->True,PlotLabel->"One Day Autocorrelation Abs[LogReturn]--Scrambled Data",Frame->True,Axes->False]

The above phenomena could occur if the time series is not stationary, but has a varying scale factor.  The γ parameter of a stable distribution is proportional to the mean absolute deviation of the data.   To explore this in greater depth we will use higher resolution data, with prices collected every minute of the trading day.  The trading days are concatenated so there is no gap between a closing price and the next open.  The SPY series here starts in July 2007.

(*data=Uncompress[Import["http://www.mathestate.com/tools/Financial/data/spy.dat","Text"]];*)
data=<<"C:\\Mathematica 7\\Portfolio\\OneMinuteData\\SPY";
lr=LogReturn[data[[All,2]]];
absecf0=Compile[{t,{s,_Real,1}},Abs[Mean[Exp[I t s]]]-1/E];
ABSECF0[t_?NumericQ,s_]:=absecf0[t,s];
SScale[s_]:=Module[{cc,c0,c1},
c1=MeanDeviation[s];
(*generally this loop should not be entered*)
While[absecf0[1/c1,s]<0,c1=2c1;Print["Upper Bound Loop"]];
c0=c1;
While[absecf0[1/c0,s]>0,c1=c0;c0=0.8 c0];(*bracket solution from above*)
cc/.FindRoot[ABSECF0[1/cc,s],{cc,  c0,c1},Method->"Brent"]
];
(*partition data by days*)
dtt=Rest[data[[All,1]]];
dtlr=Transpose[{dtt,lr}];
dtlr1=Split[dtlr,Take[#1[[1]],3]===Take[#2[[1]],3]&];
lrp=#[[All,2]]&/@dtlr1;
(*end day partition routine*)
scl=SScale/@lrp;
dtpr=Last/@Split[Rest[data][[All,{1,2}]],Take[#1[[1]],3]===Take[#2[[1]],3]&];
dt=Last[#][[1]]&/@dtlr1;
DateListPlot[dtpr,Joined->True,PlotRange->All,PlotLabel->"SPY Close"]
DateListPlot[Transpose[{dt,scl}],Joined->True,PlotRange->All,PlotLabel->"SPY Intraday Gamma",PlotStyle->Darker[Red]]

The plots above show the closing price for each day and the stable gamma calculated from the 391 minute by minute data points for each trading day.  Clearly there is a marked variation of the intraday scale factor, gamma.  Note also the opposite general slopes of the price and gamma tracings.  When we look at the autocorrelation function of the intraday data, a new level of structure is revealed.

lr1=Flatten[lrp/scl];
ListPlot[{AutoCorrelationF[Abs[lr],15 391
],AutoCorrelationF[Abs[lr1],15 391
]},PlotRange->All,Joined->True,PlotLabel->"One Day Autocorrelation Abs[LogReturn]",Frame->True,Axes->False]

The upper (blue) series is the raw data.  There are two components, one is a cyclic change that occurs at the day to day level with a sharp spike at the return caused between the sequence of  a market close and a market open where actual physical time is not continuous.  Gradually the blue curve decays over time leaving a daily autocorrelation cycle that stays close to zero after the opening spike, with a fall after the open and a rise before the next close.  The red curve is the same data but each day's data has been divided by the intraday gamma, shown in the previous plot.   This removes the day to day serial dependent structure leaving only the intraday cycle.  Although it is not easy to see this extra structure when you look at a time series of returns, it is actually easy to hear it.  The display below converts the log return series to sound using the Mathematica function ListPlay.  When you listen to the data, you will notice the crescendo of amplitude seen in November, August, and January as well as the mechanical repetitive sound related to the intraday cycle.

ListPlay[lr,SampleRate->5 391,SampleDepth->16]

CellPrint[TextCell[Hyperlink["Play Market Sound","market.wav"],"Text",TextAlignment->Center]]

MouseOver

Immediately below are the stable parameters of the minute by minute log returns and the log-log fit.  We se the same phenomenon of tails too light for the stable fit as we did in the day to day data, there are now 391 points for every day, yet the pattern is strikingly similar.

p2=StableMLFit[lr];
p3=StableMLFit[lr1];
TableForm[{p2},TableHeadings-> {None,{"α","β","γ","δ"}},TableSpacing-> {1,5}]
Quiet[DisplayLogStableFit[lr,p2]]
TableForm[{p3},TableHeadings-> {None,{"α","β","γ","δ"}},TableSpacing-> {1,5}]
Quiet[DisplayLogStableFit[lr1,p3]]

The plot and parameters immediately above shows the fit of the log returns that scaled by dividing by the daily gamma.  The division results in a gamma close to one as expected, but now the tail fit is much closer to the calculated stable fit.  The intraday cycle has been left and shows up as a obvious wiggle in the tail.  Note that the alpha parameter calculated with this data correction is now much higher.  

If we scale the intraday gamma to a daily gamma using the formula:

(Subscript[γ, daily])=(Subscript[γ, intraday]) (391^(1/α))

We find that gamma seems too high if we use the alpha from the uncorrected fit, but just about right if we use the alpha from the scaled data.

CellPrint[TextCell["Scaled by α = "<>ToString[p2[[1]]]<>", γ-daily = "<>ToString[p2[[3]] 391^(1/p2[[1]])],"Text"] ]
CellPrint[TextCell["Scaled by α = "<>ToString[p3[[1]]]<>", γ-daily = "<>ToString[p2[[3]] 391^(1/p3[[1]])],"Text"]]

Scaled by α = 1.40774, γ-daily = 0.0230707

Scaled by α = 1.75942, γ-daily = 0.00988539

At the risk of torturing the data, we can use a similar strategy to remove the intraday cycle.  Here we fit the adjusted scale factor minute by minute through the day to a parabola for the correction.  For the sample we select only full trading days; for the fit we exclude the return from the previous close to the open.

lrp1=Take[#,391]&/@Select[lrp,Length[#]≥ 391&];
minscl=SScale/@Transpose[lrp1];
y=Fit[Rest[Transpose[{Range[391],minscl}]],{x^2,x,1},x];
c=Integrate[y,{x,1,391}]/390;
q0=ListPlot[minscl/c];
q1=Plot[y/c,{x,1,391},PlotStyle->{Black,Thickness[0.008]}];
Show[q0,q1]

With this strategy we have removed all but the day to day spikes in the autocorrelation function and the oscillation in the tail fit.  The estimated parameters, {α, β} change only a little; {γ, δ} no longer have any meaning relative to our original data.

mincor=y/c/.x-> Range[391];
scl1=SScale/@lrp1;
lr2=Flatten[Transpose[Transpose[lrp1/scl1]/mincor]];
ListPlot[{AutoCorrelationF[Abs[lr2],15 391]},PlotRange->All,Joined->True,PlotLabel->"One Day Autocorrelation Adjusted Abs[LogReturn]",Frame->True,Axes->False]
p4=StableMLFit[lr2];
TableForm[{p4},TableHeadings-> {None,{"α","β","γ","δ"}},TableSpacing-> {1,5}]
Quiet[DisplayLogStableFit[lr2,p4]]

We have demonstrated some of the problems with fitting financial data to stable distributions.  Financial data are not from a stationary distribution.  There is some hope that the alpha parameter might be somewhat fixed with significantly higher value than obtained by the assumption that all the parameters are stationary, but this is not easy to prove.  A day to day variation of the scale factor confounds the problem.  It can be estimated and scaled away, but there is also a cyclic intraday pattern.  The daily variation also follows a quarterly earnings reporting cycle and is disturbed by other recurring events such as Federal Reserve Board meetings.  

Nevertheless it seems reasonable to postulate that the price formation process arises from sums of returns (or small price differences) that are driven by a fat-tailed distribution in market order books.  The apparent scale may vary by rapid changes in numbers of orders presented per unit time with a cascading effect that produces feedback in both numbers of new orders and the spacing of orders in the limit order books.  The process can quickly become complicated, but the data output by a market system suggests that sums of returns over intervals as short as a minute seem to converge to stable distributions which vary their parameters in a cyclic pattern during the trading day and with a serially dependent scale factor from day to day.  Such a system can be modelled and studied, but it does not have the convenience of simple mathematical laws that might be attributed to independent and identically distributed random variables.

Observation of the γ-parameter of the intra-day data, reveals that, although it shows a strong pattern of serial dependence, over longer time frames it might be modeled as a random variable with a lognormal distribution.  This leads to a mixture distribution, which we are studying, the LogNormalStable Distribution.   We have examined the model carefully with the most recent highly volatile data and it seems to be holding up extremely well.

Here is a link to calculate stable parameters for individual securities.

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}]],  "../r07/Volatility.html"], {5,0}]];
t3 = Graphics[Text[Hyperlink[Framed[Style["Download Notebook", {txcolor, Bold, Larger}]], "../r06/MarketData.nb"], {5, 0}]];
t4 = Graphics[Text[Hyperlink[Framed[Style["Previous", {txcolor, Bold, Larger}]], "../r05/FatTails.html"], {0, 0}]];

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

© Copyright 2010 mathestate    Wed 24 Mar 2010