Financial Time Series

This page demonstrates a Mathematica notebook for analyzing financial time series of equity prices.  The notebook may be downloaded from http://pages.suddenlink.net/rhr/fft01/FinancialTimeSeries.nb  It can be read with Mathematica or the free Wolfram CDF Player.  While everything shown in the notebook can be done directly in Mathematica, the code in the notebook uses some other packages which are not generally available.  Please contact me for the rest of the code if you would like to use the algorithms - this is an evolving work in progress.  Bob Rimmer

Load Procedures, Packages, Notebooks

NotebookEvaluate["C:\\Mathematica 9\\MarketModel2013\\S23\\S23Distribution.nb"]

ECF=Compile[{{t,_Real},{sample,_Real,1}},
(Total[Cos[ t sample]]+ I Total[Sin[ t sample]])/Length[sample],CompilationTarget->"C",Parallelization->True,RuntimeAttributes->{Listable}];

Needs["Stable`"]

(*Should also set up NS23CDF routines*)

NS23PDF[x_,{a_,b_,c_,d_}]:=6 StablePDF[x,{a/2,b,c,d}] (1-StableCDF[x,{a/2,b,c,d}]) StableCDF[x,{a/2,b,c,d}];

NS23QkPDF[x_,{a_,b_,c_,d_}]:=6 StableQkPDF[x,{a/2,b,c,d}] (1-StableQkCDF[x,{a/2,b,c,d}]) StableQkCDF[x,{a/2,b,c,d}];

tickerSymbol ="SPY";
startDate=DatePlus[Take[DateList[],3],{-10 ,"Year"}]

The time series plots on linear and logarithmic price scales.

data=FinancialData[tickerSymbol,"OHLC",startDate];
DateListPlot[Transpose[{data[[All,1]],data[[All,2]][[All,4]]}],PlotRange->All,Joined->True,PlotLabel-> tickerSymbol<>" Closing Price"]
DateListLogPlot[Transpose[{data[[All,1]],data[[All,2]][[All,4]]}],PlotRange->All,Joined->True,PlotLabel-> tickerSymbol<>" Closing Price\nLog Scale"]

Price time series often include the intraday high and low prices.  These price ranges may be converted to a standardized measure by taking the Log[High/Low].  This measure gives an estimate of the intraday volatility.

lnRange=Log[#2/#3]&@@@data[[All,2]];
DateListPlot[Transpose[{data[[All,1]],lnRange}],PlotRange->All,PlotLabel-> "Daily Log[High/Low]\n"<>tickerSymbol]
DateListPlot[Transpose[{Drop[data[[All,1]],4],MovingAverage[lnRange,5]}],PlotRange->All,Joined-> True,PlotStyle->Darker[Red],PlotLabel-> "Daily Log[High/Low]\nFive Day Moving Average "<>tickerSymbol]

These data are typically well fit by a generalized extreme value distribution.  Although volatility tends to cluster over varying time frames, a sample over a large time frame, will pick up the long term risk and tail behavior of the risk distribution.  The fit to the right tail of the 1-CDF[GEVDistribution, x] versus probability is shown below on a log-log plot where the tail behavior of the Log[Range] typically becomes linear at the extreme.  The GEV distribution in Mathematica is called the MaxStableDistribution[μ, σ, ξ] -- the parameters are shown below.  In the plot the fit distribution is in red and the empirical distribution is in blue, showing the closeness of the fit.  Aside from providing a useful method to model volatility over time, the statistical analysis of this measure as an intraday extreme may shed some light on the power tail nature of the underlying price formation process.  If the extremes measured by the log ratio were being generated by a normal distribution process or even heavier tailed process like an exponential distribution, the parameter, ξ, should be close to zero.  Instead ξ is almost always securely in the Frechet domain of attraction.  The value is usually less than 1/2 suggesting that the distribution on summation should eventually lead to a normal distribution, however the underlying process still has far heavier tails than an exponential distribution.

rngD=EmpiricalDistribution[lnRange];
gev=EstimatedDistribution[lnRange,MaxStableDistribution[μ,σ,ξ]]
LogLogPlot[{1-CDF[rngD,x],1-CDF[gev,x]},{x,0.0025,Max[lnRange]},PerformanceGoal->"Speed",PlotRange->All,Frame->True,GridLines->Automatic,FrameLabel->{"Log[High/Low]","1 - Probability"},PlotLabel->"Log[high/low] GEV Fit\nGEV Tail α: "<>ToString[1/gev[[3]],TraditionalForm]]

Extreme Event Clustering

The graph below shows similar information in a different format.  The dots are sequential daily logarithmic returns rather than the intraday measure Log[High/Low].  The dots in red are outside the probability range (0.02, 0.98).  The main idea of this graph is to clearly show that extreme events are not evenly distributed across time as would be expected if the underlying price formation process were random.  Instead there is likely a feedback loop in the process so that extremes tend to cluster.  The point is that although it is possible to fit these data collected over a long period of time rather well to a particular distribution, the probabilities of events from that distribution will not have the same predictive value as they would if the events were being generated by a random process.  So what follows should be used in practice with extreme caution.  If you are operating in an epoch when extreme events are not occurring the procedures will over estimate risk.  If you are in an epoch during which extreme events are occurring, risk will still significantly be underestimated because of the clustering of extreme events where an ordinary random process would expect extremes to be randomly distributed.

lr=Differences[Log[data[[All,2]][[All,4]]]];
dt=Rest[data[[All,1]]];
pLevel=0.02;
lrQrng=Quantile[lr,{pLevel,1-pLevel}];
DateListPlot[Tooltip@{Transpose[{dt,lr}],Select[Transpose[{dt,lr}],!lrQrng[[1]]<#[[2]]<lrQrng[[2]]&]},PlotStyle-> {Darker[Blue],Red},PlotRange->All,PlotLabel-> "Daily Log Returns\np = "<>ToString[pLevel,TraditionalForm]<>" Extremes Red"]

Stable 2:3 Order Distribution Model Fit

This distribution is described in more detail here: http://pages.suddenlink.net/rhr/fin05/MarketThoughtExperiment.html.  This distribution tends to fit the tail behavior of stock market logarithmic returns very well.  The distribution tail exponent, α, has support in a range (0, 4), but at exactly 4 it picks up the tail behavior of a normal distribution.  Usually for U.S. market returns, α is in the 2.5 to 3.5 range.  With a tail exponent at this level, sums of returns should eventually yield a normal distribution by the central limit theorem, but the convergence is slow, and in general convergence will not occur over even long time frames because of the non random clustering of extreme events.  Financial returns are neither identically nor independently distributed over time.  And this slow convergence aspect of market tail behavior by the central limit theorem is lost entirely if one mistakenly assumes an underlying normal distribution.  It is possible to use the FFT to calculate the convergence over varying intervals, but doing this requires that the assumed market behavior remain similar to the data epoch and that events occur independently.  The next section of this notebook will demonstrate how to do this, but the main value in the exercise is to show the effect of the tail behavior in a randomized setting rather than in the real world.

The empirical distribution will be compared to the fit to the stable 2:3 order statistic distribution.  This is the distribution which arises from repeatedly selecting the median of three stable random variables.  The rationale underlying the use of such a distribution is that it can be seen in the order books of the continuous double auction that limit bid and ask prices contain frequent jumps which show heavy power tails.  A stable distribution can conveniently model very heavy tailed behavior but if it is directly applied to logarithmic returns from financial time series, it picks tails which are heavier than the data.  The idea here is that many of the orders placed in the price formation process in the continuous double auction never result in a transaction.  But what we see in reported prices result from transactions between buyers and sellers.  The median order statistic distribution from a stable distribution with maximum rank three provides way to select a subset of random variables representing a compromise between buyers and sellers.  It may not be a direct process in price formation, but may ultimately result from many transactions from heavy Pareto tails in the order books summing to something like a stable distribution when looked at over time frames longer than a minute.  The selection of this particular subset of stable random variables has the property that it allows distribution tail behavior which has precisely twice the magnitude of the tail behavior of the underlying stable distribution.  It is also easy to generate these random variables for simulation.

lrD=EmpiricalDistribution[lr];

parm0=Quiet@S23MLFit[lr,{3,0,0.01,0}];//AbsoluteTiming
parm0

The parameters above are in the format {α, β, γ, δ}, where α is the actual tail behavior of the distribution function and the other parameters correspond to the parent stable distribution which has an α which is half the value shown above.  The plot below shows the right tail in red and the left tail in blue.  The jagged curves are the empirical distribution.  In this case the actual empirical data tail is a little lighter (higher magnitude α, or steeper negative slope on the plot) than that of the smoother fit curves.

LogLogPlot[{CDF[lrD,-x],1-CDF[lrD,x],CDF[S23OrderDistribution[Sequence@@parm0],-x],1-CDF[S23OrderDistribution[Sequence@@parm0],x]},{x,0.001,0.5},PlotRange->All,PerformanceGoal->"Speed",Frame->True,FrameLabel->{"Log[High/Low]","Probability"},GridLines -> Automatic,PlotStyle->{Darker[Blue],Darker[Red],Darker[Blue],Darker[Red]},PlotLabel-> "Tail Fit with Empirical Distribution"]

FFT Model and Empirical Convolution

FFT Vectors and Functions

(*the procedures which follow are set up to use the default parameters for the Mathematica FFT functions*)
(* the rationale for this choice rather than the parameters for statistics is to be able to use the code also in the CUDAFourier function*)

XVector  = Function[{SR, nn},Module[{n},
n=2 Floor[nn/2];
Range[-n/2,n/2-1,1]/(SR)]];

ReOrder=With[{n=Length[#]},Join[Drop[#,Floor[n/2]],Drop[#,-Floor[n/2]]]]&;

TVector=Function[{SR, nn ,tlen},Module[{n0,n1},
n0=2 Floor[nn/2];n1=2 Floor[tlen/2];
2 Pi  SR /n0 Range[0.,n1]]];

(*SR is the sampling rate or the number of sample intervals to be taken per unit of the scaleRange*)
(*nn is the total sample size.  For FFT having this a power of 2 should improve efficiency.*)

nn=2^15
SR=2^10
scaleRange=nn/SR

xvec=XVector[SR,nn];
xvec[[1]]//N
xvec[[-1]]//N

(*For stable and empirical characteristic functions the value of the characteristic function at minus t is the conjugate of the value at t, so only half the domain of t carries information.*)

(*{alpha,beta}={1.09,0.8};*)
n1=2^14;
tvec=TVector[SR,nn,n1];

xvec[[nn/2+n1]]//N

(*For summing random variables it makes sense to remove the trend by subtracting  the mean*)

mlr=Mean[lr]

cfvec=ECF[tvec,lr-mlr];//AbsoluteTiming

Length[cfvec]

cf1=ReOrder[Join[ConstantArray[0,nn/2-n1],Conjugate[Reverse[Drop[cfvec,1]]],Drop[cfvec,-1],ConstantArray[0,nn/2-n1]]];

Length[cf1]==nn

cf1[[1]]

(*Construct distribution of sum of 21 random variables from cf1*)

alternateSign=Riffle[ConstantArray[1,nn/2],ConstantArray[-1,nn/2]];

(* SR/Sqrt[nn ] scales for default Fourier settings*)

pdf21=    Re[InverseFourier[  (SR/Sqrt[nn ]) cf1^21 alternateSign  ]];

(*The FFT tries to make the output function begin and end at the same point.  But there is also some curvature in the tails at the very end when the distribution is symmetric.  Most likely this is due in part to filtering of the high frequency harmonics.*)

(*one way to have the end and beginning points closer might be to recenter the distribution, using the 0-parameterization.*)

This next section uses the empirical characteristic function (ecf) derived directly from the data points, where the ecf is defined below and parameter, n, represents the sample size and (Subscript[X, k]) is an individual random variable.  i = (Sqrt[-1]).

ecf(t)=(1/n) (∑_(k=1)^n)(e^(i  t (Subscript[X, k])))

The ecf is sampled at equally spaced points over a large domain and the inverse FFT can convert this to a sampled density function for the data distribution.  An explicit formula for the stable 2:3 order statistic distribution is not known, by me at least, and may not exist, but it can be approximated by sampling the density and using the FFT to generate a sampled characteristic function.  Both of these sampled characteristic functions can be scaled to a distribution which is the sum of some number, m, random variables by raising each element of the sampled set to the m th power.  In the example m will be 21 to represent 21 trading days in a month, converting the daily distribution of time series of logarithmic returns to a monthly time series.  After raising each characteristic function set to the power, it can be inverted back to a sampled density function from which probability calculations may be made, remembering the caveats about the degree of Independence and stationarity of the data.

Code

mx=Max[pdf21]
pos=Position[pdf21,mx][[1]][[1]]

xvec[[pos]]

ListLogLogPlot[{{-#1,#2}&@@@Take[Transpose[{xvec,pdf21}],pos],Take[Transpose[{xvec,pdf21}],{pos+1,nn}]},GridLines-> Automatic,PlotRange->{{0.01,5},All}]

(*So it takes very small amount of information to construct most of the pdf*)

(*Need an algorithm to select where to take the end point.*)

smoothLimit=1.;
ipdf21=Interpolation[Select[Transpose[{xvec,pdf21}],Abs[#[[1]]]<smoothLimit&](*,Method->"Spline"*)]

epts={ipdf21[[1]][[1]][[1]],ipdf21[[1]][[1]][[-1]]}

NIntegrate[ipdf21[x],{x,epts[[1]],epts[[2]]}]

(*Stable 2:3 order distribution fit mean adjusted*)

parm=Quiet@S23MLFit[lr-mlr,{3,0,0.01,0}]

pdf1=    Re[InverseFourier[  (SR/Sqrt[nn ]) cf1 alternateSign  ]];

ipdf1=Interpolation[Select[Transpose[{xvec,pdf1}],Abs[#[[1]]]<0.3&](*,Method->"Spline"*)]

The density curve derived from the FFT is shown below in red and compared to the density of the fit are shown below.  In this case mean adjusted data are used to remove any trend in the data from the convolution process.  There are troublesome oscillations in the FFT density curve and because of these, it is not positive definite as would be expected for a density function, but interestingly much of this artifact disappears on convolution.

(*These are slower much slower Mathematica based procedures*)(*Plot[{PDF[S23OrderDistribution[parm],x],NS23PDF[x,parm]},{x,-0.1,0.1},PlotRange->All]
Plot[{PDF[S23OrderDistribution[parm],x]-NS23PDF[x,parm]},{x,-0.1,0.1},PlotRange->All]*)
Plot[{NS23PDF[x,parm],ipdf1[x]},{x,-0.1,0.1},PlotRange->All]

PDF Tail behavior.

This compares model tail to data tail.  There basically are not enough points to tell whether the model is good.  The tails of the density function are shown on the log-log plot, using the absolute value of negative x.  The blue and red are the left and right tails of the fit and the yellow and green are the left and right tails of the FFT derived density.

LogLogPlot[{NS23PDF[-x,parm],NS23PDF[x,parm],ipdf1[-x],ipdf1[x]},{x,0.001,0.3},PlotRange->All,PerformanceGoal->"Speed"]

Scale the S23 to 21 days.

Length[cfvec]

n1

Length@Take[xvec,{nn/2-n1,nn/2+n1}]

(*The Nolan Qk routines work well here and very fast*)

pdfvec=NS23QkPDF[xvec,parm];//AbsoluteTiming

Need to test the accuracy of the Mathematica built in grid, but there is no CDF grid to make s23 distribution, Could test the grids in StableM7`, in the past these were set to be a little more accurate for alpha > 1.

s23cf=Fourier[pdfvec];//AbsoluteTiming

Length[s23cf]

s23cf[[1]]

nn

Length[s23cf]

ECF and Model CF Comparison

The plot below shows the points of the ecf and the FFT derived characteristic function points projected onto the complex plane with the real values on the x-axis and the imaginary values on the y-axis.  The FFT generated characteristic function is free of the oscillations which are present in the empirical distribution.  The range of a statistical characteristic function is confined to the unit circle in the complex plane.

ListPlot[{{Re[#],Im[#]}&/@cf1,{Re[#],Im[#]}&/@(alternateSign s23cf/s23cf[[1]])},PlotRange->All,Frame->True,FrameLabel->{"Re","Im"},PlotLabel->"Empirical CF Blue\ns23 model FFT CF Red"]

After raising every point in the complex plane above to the 21st power and inverting the FFT, an interpolated density curve is plotted below.  The blue and red curves respectively belong to the ecf and FFT calculated characteristic functions.

(*alternateSign added in case of even numbered powers*)
(*interestingly this doesn't need any imaginary numbers removed*)
s23pdf21=InverseFourier[(SR/Sqrt[nn ]) alternateSign (alternateSign s23cf/s23cf[[1]])^21];
is23pdf21=Interpolation[Transpose[{xvec,s23pdf21}]];
Plot[{ipdf21[x],is23pdf21[x]},{x,-0.29,0.29}]

On the log-log plot of the four density tails the colors are changed.  The red and blue belong to the FFT red curve above and the yellow and green tails belong to the blue curve ecf data.  The power tail behavior of the empirical distribution has been lost on convolution because the empirical data did not have enough data points on in the tails of the actual data sample to preserve it.  The stable 2:3 order statistic distribution was used with enough sampled points to preserve the linearity on a log-log plot of the power tail distribution.  As convolution is performed with greater number of convolutions, the linear portion of the curve occurs ever further out on the tails until it disappears after infinitely many convolutions.  The model can be used to determine at what point in the distribution the power tails break away.

LogLogPlot[{is23pdf21[-x],is23pdf21[x],ipdf21[-x],ipdf21[x]},{x,0.01,2},PlotRange->All,PerformanceGoal->"Speed"]

There is artifact which occurs with the FFT and that is due to aliasing when the two end points of the sampled function are not the same.  To avoid having this be problematic a very large sample is used for the FFT, so the aliasing wrap around effect is not significant.  The left tail is blue and the right red.  The values on the abscissa of both this plot and the plot above are the absolute values of x and the ordinate axis shows the density values.

(*The very distal tails demonstrate the Fourier wrap around*)
ListLogLogPlot[{{-#1,#2}&@@@Take[Transpose[{xvec,s23pdf21}],2000],Take[Transpose[{xvec,s23pdf21}],-2000]},PlotLabel->"Extreme FFT Density Points"]

Tail fit: right tail

Fit[Log[Select[Transpose[{xvec,s23pdf21}],15.2<#[[1]]<15.6&]],{x,1},x]

Tail fit : left tail

Fit[Log[{-#1,#2}&@@@Select[Transpose[{xvec,s23pdf21}],-15.6<#[[1]]<-15.2&]],{x,1},x]

Expected PDF tail exponent from ML fit

parm[[1]]+ 1

Finally below is a price probability curve using the current price and the one month distribution of logarithmic returns scaled from daily price returns using the FFT of the model density -> characteristic function ^ 21 -> inverse FFT density -> integrated to a distribution function.

pr=FinancialData[tickerSymbol];
Off[NIntegrate::nlim];
Plot[NIntegrate[is23pdf21[x],{x,xvec[[1]],Log[t/pr]}],{t,150,230},GridLines->Automatic,Frame->True,PerformanceGoal->"Speed",PlotLabel-> "One Month Price Probability Curve\nCurrent Price "<>ToString[pr,TraditionalForm]]
On[NIntegrate::nlim];

Financial Data Analysis Home

CellPrint[TextCell["\n\n© Copyright " <> ToString[DateList[][[1]]] <> " Robert H. Rimmer, Jr.    " <> DateString[DateList[], "DateShort"], "Text"]]

© Copyright 2014 Robert H. Rimmer, Jr.    Wed 28 May 2014