Convolution

Needs["StableM`"];
Needs["Histograms`"];
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]))]]];

Convolution in the transform domain is simple multiplication.  Here is an example of convolution of symmetric stable distribution with a square wave, the amplitude is adjusted so that the area under the curve is 1, i.e. it is a statistical density function.

square[x_, w_] := (Sign[x + w/2] - Sign[x - w/2])/(2 w);
Plot[square[x, 1], {x, -2, 2}]

square[x_,w_]:= (Sign[x+w/2]-Sign[x-w/2])/(2 w);
Plot[square[x,1],{x,-2,2}]

FourierTransform[square[x, w], x, t, Assumptions -> w ∈ Reals]
Limit[-1/(t w) Sqrt[2/Pi] Sin[t w/2] (-1 + UnitStep[-w] UnitStep[w]), t -> 0]

FourierTransform[square[x,w],x,t,Assumptions-> w∈ Reals]
Limit[-1/(t w) Sqrt[2/Pi] Sin[t w/2] (-1+UnitStep[-w] UnitStep[w]),t->0]

squaref[t_, w_] := -1/(t w) Sqrt[2/Pi] Sin[t w/2] (-1 + UnitStep[-w] UnitStep[w]);
squaref[0, w_] := -(-1 + UnitStep[-w] UnitStep[w])/Sqrt[2 Pi];
scale = 8192/16;
sqfilt =  scale  DFTCF[squaref[#1, 5] &, 16, 8192];
ListPlot[Transpose[{XTable[16, 8192], ReOrder[sqfilt]}], PlotRange -> {{-scale/2, scale/2}, All}]

squaref[t_,w_]:=-1/(t w) Sqrt[2/Pi] Sin[t w/2] (-1+UnitStep[-w] UnitStep[w]);
squaref[0,w_]:=-(-1+UnitStep[-w] UnitStep[w])/Sqrt[2 Pi];
scale=8192/16;
sqfilt=  scale  DFTCF[squaref[#1,5]&,16,8192];
ListPlot[Transpose[{XTable[16,8192],ReOrder[sqfilt]}],PlotRange->{{-scale/2,scale/2},All}]

scf = DFTCF[CF1c[#1, 1.3`, 0, 1, 0] &, 16, 8192];
ListPlot[Transpose[{XTable[16, 8192], scale ReOrder[scf]}], PlotRange -> {{-scale/2, scale/2}, All}]

scf= DFTCF[CF1c[#1,1.3`,0,1,0]&,16,8192];
ListPlot[Transpose[{XTable[16,8192],scale ReOrder[scf]}],PlotRange->{{-scale/2,scale/2},All}]

conv = Chop[InverseFourier[ scf sqfilt]];
g1 = ListPlot[Transpose[{XTable[16, 8192], ReOrder[conv]}], PlotRange -> {{-20, 20}, All}, Joined -> True, PlotStyle -> Blue];
conv = Chop[InverseFourier[ scf]];
g2 = ListPlot[Transpose[{XTable[16, 8192], ReOrder[conv]}], PlotRange -> {{-20, 20}, All}, Joined -> True];
conv = Chop[InverseFourier[sqfilt/scale]];
g3 = ListPlot[Transpose[{XTable[16, 8192], ReOrder[conv]}], PlotRange -> {{-20, 20}, All}, Joined -> True];
Show[g2, g3, g1]

conv=Chop[InverseFourier[ scf sqfilt]];
g1=ListPlot[Transpose[{XTable[16,8192],ReOrder[conv]}],PlotRange->{{-20,20},All},Joined->True,PlotStyle->Blue];
conv=Chop[InverseFourier[ scf]];
g2=ListPlot[Transpose[{XTable[16,8192],ReOrder[conv]}],PlotRange->{{-20,20},All},Joined->True];
conv=Chop[InverseFourier[sqfilt/scale]];
g3=ListPlot[Transpose[{XTable[16,8192],ReOrder[conv]}],PlotRange->{{-20,20},All},Joined->True];
Show[g2,g3,g1]

Convolution of two characteristic functions will give the characteristic function of the the sum of a random variable from each function.  The following will give the density of the sum of an S(1.3, 1, 1, 0;1) and a S(2 ,0, 1, 0;1) random variable.

scale = 8192/16;
skew = DFTCF[CF1c[#, 1.3, 1, 1, 0] &, 16, 8192];
nrml = scale DFTCF[CF1c[#, 2, 0, 1, 0] &, 16, 8192];
conv1 = InverseFourier[ skew nrml] // Chop;
g0 = ListPlot[Transpose[{XTable[16, 8192], ReOrder[conv1]}], PlotRange -> {{-50, 50}, All}, Joined -> True, PlotStyle -> Thickness[0.008]]

scale=8192/16;
skew=DFTCF[CF1c[#,1.3,1,1,0]&,16,8192];
nrml=scale DFTCF[CF1c[#,2,0,1,0]&,16,8192];
conv1=InverseFourier[ skew nrml]//Chop;
g0=ListPlot[Transpose[{XTable[16,8192],ReOrder[conv1]}],PlotRange->{{-50,50},All},Joined->True,PlotStyle->Thickness[0.008]]

The convolution can be obtained very slowly with numerical integration.

AbsoluteTiming[g1 = Plot[NIntegrate[SPDF[t, {1.3, 1}, 1] SPDF[x - t, {2, 0}, 1], {t, -∞, ∞}, Method -> {Automatic, "SymbolicProcessing" -> 0}], {x, -50, 50}, PlotStyle -> Orange, PlotRange -> All];]
Show[g0, g1]

AbsoluteTiming[g1=Plot[NIntegrate[SPDF[t,{1.3,1},1] SPDF[x-t,{2,0},1],{t,-∞,∞},Method-> {Automatic,"SymbolicProcessing"->0}],{x,-50,50},PlotStyle->Orange,PlotRange->All];]
Show[g0,g1]

Obtain the result by simulation.

s = SRandom[10000, {1.3, 1}] + SRandom[10000, {2, 0}];
g2 = Histogram[s, HistogramScale -> 1, HistogramRange -> {-10, 30}];
Show[g2, g0, PlotRange -> {{-10, 30}, All}]

s=SRandom[10000,{1.3,1}]+SRandom[10000,{2,0}];
g2=Histogram[s,HistogramScale->1,HistogramRange->{-10,30}];
Show[g2,g0,PlotRange->{{-10,30},All}]

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}]], "../ft06/EmpiricalCF.html"], {0, -0.5}]];
t3 = Graphics[Text[Hyperlink[Framed[Style["Download Notebook", {txcolor, Bold, Larger}]], "../ft05/Convolution.nb"], {5, 0}]];
t4 = Graphics[Text[Hyperlink[Framed[Style["Previous", {txcolor, Bold, Larger}]], "../ft04/Differentiation.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 29 Feb 2008