MethodolComputApplProbabmanuscriptNo. (willbeinsertedbytheeditor) An Analytic Expression for the Distribution of the Generalized Shiryaev–Roberts Diffusion TheFourierSpectralExpansionApproach AlekseyS.Polunchenko · GrigorySokolov Received:date/Accepted:date 6 1 0 Abstract Weconsiderthequickestchange-pointdetectionproblemwheretheaimistode- 2 tect the onset of a pre-specified drift in “live”-monitored standard Brownian motion; the n change-pointisassumedunknown(nonrandom).Thetopicofinterestisthedistributionof a the Generalized Shryaev–Roberts (GSR) detection statistic set up to “sense” the presence J ofthedrift.Specifically,wederiveaclosed-formformulaforthetransitionprobabilityden- 5 sity function (pdf) of the time-homogeneous Markov diffusion process generated by the 1 GSRstatisticwhentheBrownianmotionundersurveillanceis“drift-free”,i.e.,inthepre- change regime; the GSR statistic’s (deterministic) nonnegative headstart is assumed arbi- ] E trarilygiven.Thetransitionpdfformulaisfoundanalytically,throughdirectsolutionofthe M respectiveKolmogorovforwardequationviatheFourierspectralmethodtoachievesepara- tionofthespacialandtemporalvariables.Theobtainedresultgeneralizesthewell-known . t formulaforthe(pre-change)stationarydistributionoftheGSRstatistic:thelatter’sstation- a t ary distribution is the temporal limit of the distribution sought in this work. To conclude, s weexploittheobtainedformulanumericallyandbrieflystudythepre-changebehaviorof [ the GSR statistic versus three factors:(a) drift-shift magnitude, (b) time, and (c) the GSR 1 statistic’sheadstart. v 8 Keywords GeneralizedShiryaev–Robertsprocedure·Kolmogorovforwardequation· 6 Markovdiffusionprocesses·Quickestchange-pointdetection·Sequentialanalysis 8 3 0 A.S.Polunchenko 1. DepartmentofMathematicalSciences StateUniversityofNewYorkatBinghamton 0 Binghamton,NY13902–6000,USA 6 Tel.:+1-607-777-6906 1 Fax:+1-607-777-2450 : v E-mail:[email protected] i G.Sokolov X DepartmentofMathematicalSciences r StateUniversityofNewYorkatBinghamton a Binghamton,NY13902–6000,USA Tel.:+1-607-777-4239 Fax:+1-607-777-2450 E-mail:[email protected] 2 PolunchenkoandSokolov MathematicsSubjectClassification(2000) MSC62L10·MSC60G10·MSC62M15· MSC60J60 1 Introduction Sequential(quickest)change-pointdetectionisconcernedwiththedevelopmentandevalu- ationofreliablestatisticalproceduresforearlydetectionofunanticipatedchangesthatmay (or may not) occur online in the characteristics of a “live”-monitored (random) process. Specifically, the latter is observed continuously with the intent to “flag an alarm” in the event (and as soon as) the behavior of the process starts to suggest the process may have (been)statisticallychanged.Thealarmistobeflaggedasquicklyasispossiblewithinaset tolerablelevelofthe“falsepositive”risk.See,e.g.,Shiryaev(1978),BassevilleandNiki- forov(1993),PoorandHadjiliadis(2009),VeeravalliandBanerjee(2013),Tartakovskyetal (2014)andthereferencestherein. Achange-pointdetectionprocedureisidentifiedwithastoppingtime,T,thatisadapted tothefiltration,(Ft)t(cid:62)0,generatedbytheobservedprocess,(Xt)t(cid:62)0;thesemanticsofT is thatitconstitutesaruletostopanddeclarethatthestatisticalprofileoftheobservedprocess mayhave(been)changed.A“good”(i.e.,optimalornearlyoptimal)detectionprocedureis onethatminimizes(ornearlyminimizes)thedesireddetectiondelaypenalty-function,sub- jecttoaconstraintonthefalsealarmrisk.Foranoverviewofthemajoroptimalitycriteria, see,e.g.,TartakovskyandMoustakides(2010),PolunchenkoandTartakovsky(2012),Pol- unchenkoetal(2013),VeeravalliandBanerjee(2013)or(Tartakovskyetal2014,PartII). This work focuses on the popular minimax setup of the basic change-point detection problemwheretheobservedprocess,(Xt)t(cid:62)0,isstandardBrownianmotionthatatanun- known (nonrandom) time moment ν—referred to as the change-point—may (or may not) experienceanabruptandpermanentchangeinthedrift,fromavalueofzeroinitially,i.e., E[dX ] = 0 for t ∈ [0,ν], to a known value µ (cid:54)= 0 following the change-point, i.e., t E[dX ] = µ for t ∈ (ν,∞). This is illustrated in Figure 1. The goal is to find out—as t quicklyasispossiblewithinanapriorisetlevelofthe“falsepositive”risk—whetherthe driftoftheprocessisnolongerzero.See,e.g.,PollakandSiegmund(1985),Shiryaev(1996, 2002),Moustakides(2004),FeinbergandShiryaev(2006),andBurnaevetal(2009). Moreformally,undertheaboveBrownianmotionchange-pointscenario,theobserved process,(Xt)t(cid:62)0,isgovernedbythestochasticdifferentialequation(SDE): dX =µ1l dt+dB , t(cid:62)0, withX =0, (1) t {t>ν} t 0 where(Bt)t(cid:62)0isstandardBrownianmotion(i.e.,E[dBt]=0,E[(dBt)2]=dt,andB0 = 0),µ(cid:54)=0istheknownpost-changedriftvalue,andν ∈[0,∞]istheunknown(nonrandom) change-point;hereandonward,thenotationν =0(ν =∞)istobeunderstoodasthecase whenthedriftisineffectabinitio(ornever,respectively). Toperformchange-pointdetectionundermodel(1),thestandardapproachhasbeento employ Page’s (1954) Cumulative Sum (CUSUM) “inspection scheme”. This choice may bejustifiedbythefact(establishedbyBeibel1996,byShiryaev1996andbyMoustakides 2004)thattheCUSUMschemeisstrictlyminimax-optimalinthesenseofLorden(1971); the discrete-time equivalent of this result was first established by Moustakides (1986), al- thoughanalternativeproofwasalsolaterofferedbyRitov(1990)whousedagame-theoretic argument. However, when one is interested in minimax optimality as defined by Pollak (1985), a sensible alternative to using the CUSUM scheme would be to devise the Generalized ExactDistributionoftheGeneralizedShiryaev–RobertsDiffusion 3 Standard Brownian Motion + Drift Function Drift Function a t a D Change-Point(Unknown) Standard Brownian Motion (No Drift) Time Fig.1:StandardBrownianmotiongainingapersistentdrift. Shiryaev–Roberts(GSR)procedure.ThelatterisduetoMoustakidesetal(2011)andisa headstarted version of the classical quasi-Bayesian Shiryaev–Roberts (SR) procedure that emergedfromtheindependentworkofShiryaev(1961,1963)andthatofRoberts(1966). WithPollak’s(1985)definitionofminimaxoptimalityinmind,themotivationtopreferthe GSR procedure over the CUSUM chart stems from the results obtained (for the discrete- timeanalogueoftheproblem)byTartakovskyandPolunchenko(2010);Polunchenkoand Tartakovsky (2010) and then also by Tartakovsky et al (2012) who showed that the GSR procedure with a carefully designed headstart may be faster (in Pollak’s 1985 sense) than theCUSUMscheme;asamatteroffact,TartakovskyandPolunchenko(2010);Polunchenko andTartakovsky(2010)provedtheGSRprocedure(witha“finetuned”headstart)tobethe fastest(i.e.,thebestonecando)intwospecific(discrete-time)scenarios.Foranattemptto extendtheseresultstotheBrownianmotionscenario,see,e.g.,Burnaev(2009). Morespecifically,theGSRprocedurecallsforstoppingassoonastheGSRdetection statistic,(Rtr)t(cid:62)0,hitsacriticallevelknownasthedetectionthreshold.Thelatterissetso ashavethe“falsepositive”riskatadesired“height”.LetP (P )denotetheprobability ∞ 0 measure(distributionlaw)generatedbytheobservedprocess,(Xt)t(cid:62)0,undertheassump- tionthatν = ∞(ν = 0);notethatP istheWienermeasure.Let P | (P | )bethe restrictionofprobabilitymeasureP ∞(P )tothefiltrationF .Further,∞deFfitne 0 Ft ∞ 0 t dP | Λ (cid:44) 0 Ft , t(cid:62)0, t dP | ∞ Ft i.e.,theRadon–Nikodýmderivativeof P | withrespectto P | .Itiswell-knownthat 0 Ft ∞ Ft fortheBrownianmotionscenariounderconsideration (cid:26) µ2 (cid:27) Λ =exp µX − t , sothatdΛ =µΛ dX , Λ =1; t t t t t 0 2 cf. Shiryaev (1999) and Liptser and Shiryaev (2001). The process {Λt}t(cid:62)0 manifests the likelihood ratio to test the hypothesis H : ν = 0 against the alternative H : ν = ∞, 0 ∞ 4 PolunchenkoandSokolov andisthekeyingredientoftheCUSUMstatisticaswellasoftheGSRstatistic,(Rtr)t(cid:62)0. Specifically,tailoredtotheBrownianmotionscenarioathand,theGSRstatistic,(Rtr)t(cid:62)0, isoftheform (cid:90) t Λ Rr (cid:44)rΛ + t ds t t Λ 0 s (2) (cid:26) µ2t(cid:27) (cid:90) t (cid:26) µ2(t−s)(cid:27) =rexp µX − + exp µ(X −X )− ds, t(cid:62)0, t t s 2 2 0 where Rr = r (cid:62) 0 is the headstart (a deterministic point selected so as to optimize the 0 GSRprocedure’sperformance;see,e.g.,TartakovskyandPolunchenko2010;Polunchenko and Tartakovsky 2010; Moustakides et al 2011; Tartakovsky et al 2012; Polunchenko and Sokolov2014).WhenRr = r = 0,itissaidthatthereisnoheadstart.TheGSRstatistic, 0 (Rtr)t(cid:62)0,withnoheadstartisequivalenttotheclassicalSRstatistic.Consequently,theGSR procedurewhosestatistichasnoheadstartisequivalenttotheclassicalSRprocedure.Hence, thelabels“GeneralizedSRstatistic”and“GeneralizedSRprocedure”,whichappeartohave bothbeencoinedbyTartakovskyetal(2012). Wearenowinapositiontoformulatethespecificproblemaddressedinthispaper:to obtain an explicit closed-form formula for p (y,t|x,s) (cid:44) dP (Rr (cid:54) y|Rr = x)/dy, ∞ ∞ t s 0 (cid:54) s < t < ∞, x,y (cid:62) 0, i.e., for the P -transition probability density function (pdf) ∞ oftheGSRstatistic(Rtr)t(cid:62)0 withtheheadstart,R0r = r (cid:62) 0,assumedgiven.Thatis,in this paper the GSR statistic, (Rtr)t(cid:62)0, is effectively let “run loose” over the entire space (Rr,t) ∈ [0,∞)×[0,∞) with no detection threshold imposed, and the goal is to find t the exact transition pdf of the one-dimensional Markov diffusion process (Rtr)t(cid:62)0 under probabilitymeasureP .Morespecifically,observethat,byItôformula,theP -differential ∞ ∞ of the GSR diffusion (Rtr)t(cid:62)0 is dRtr = dt + µRtrdBt, where R0r = r (cid:62) 0, whence (Rtr)t(cid:62)0 is seen to be time-homogeneous. Therefore, p∞(y,t|x,s) depends on s and t onlythroughthedifferencet−s (cid:62) 0,anditsufficestofindp (x,t|r) (cid:44) p (x,t|r,0), ∞ ∞ x,r,t(cid:62)0.Hence,weshallconcentrateonfindingp (x,t|r),x,r,t(cid:62)0.Also,tolighten ∞ thenotation,fromnowonweshalldropthe“∞”inthesubscriptofp (x,t|r)andsimply ∞ writep(x,t|r). Theprincipalapproachweintendtoundertaketofindp(x,t|r)consistsinsolvingdi- rectlytherespectiveKolmogorovforwardequation ∂ p(x,t|r)=− ∂ p(x,t|r)+ µ2 ∂2 (cid:2)x2p(x,t|r)(cid:3), x,t,r(cid:62)0, (3) ∂t ∂x 2 ∂x2 subjectto(a)the(natural)normalizationconstraint (cid:90) ∞ p(x,t|r)dx=1 (4) 0 validforallr,t(cid:62)0,and(b)theinitialconditionlim p(x,t|r)=δ(x−r)validforall t→0+ x; here and onward δ(z) denotes the Dirac delta function so that “lim p(x,t|r) = t→0+ δ(x − r)” is to be understood as equality of distributions. We note also that naturally p(x,t|r) (cid:62) 0forallr,x,t (cid:62) 0.Theinitialconditionandthenormalizationconstraint(4) arealsotobecomplementedbytwoboundaryconditionsofspacial(i.e.,inthexvariable) type.TheseconditionsareprovidedinSection3whichisalsothesectionwherewesolve(3) explicitlyandthusobtainthemainresultofthispaper. Morespecifically,torecoverp(x,t|r)from(3)wedevisetheso-calledFourierspectral method, a separation-of-variables-type technique the general theory for which was devel- oped by Weyl (1910) and by Titchmarsh (1962). See also, e.g., Levitan (1950), McKean ExactDistributionoftheGeneralizedShiryaev–RobertsDiffusion 5 (1956),DunfordandSchwartz(1963),LevitanandSargsjan(1975).TheWeyl–Titchmarsh theorywasthenappliedandsharpenedfurtherinthecontextsofprobabilityandstochastic processes by Kac (1951, 1959), Feller (1952), Itô and McKean (1974). With the specifics deferred to Section 3, the Fourier method allows to separate the spacial variable, x, and the temporal variable, t, and thus reduce the original equation (3) to a (singular) Sturm– Liouvilleproblem.Thefundamentalsolutionsofthelatter,inturn,allowtoexplicitlycon- structtheGreen’sfunctionassociatedwiththecorrespondingSturm–Liouville(linear,dif- ferential)operator.AstheresolventoftheSturm–Liouvilleoperator,theGreen’sfunction providesanexhaustivecharacterizationoftheoperator’sproperties,includingitsspectrum. Inparticular,theGreen’sfunctionisdirectlyconnectedtothetransitionpdfp(x,t|r):the formeristheLaplacetransform(takenwithrespecttotime)ofthelatter.Therefore,getting p(x,t|r)isamatterofinvertingits(temporal)LaplacetransformgivenbytheGreen’sfunc- tion.TheinversioncanaccomplishedbyvirtueofCauchy’sResidueTheorem.Thepoles(or branchcuts)oftheGreen’sfunctionyieldtheeigenvaluesoftheSturm–Liouvilleoperator andthecorrespondingresiduesdeterminethecontributionsoftheeigenfunctionstothepdf p(x,t|r). Theformulaforp(x,t|r)foundusingtheFouriermethodinSection3isageneralization oftheresultofShiryaev(1961,1963)whoobtainedthelimitρ(x) (cid:44) lim p(x,t|r), t→+∞ i.e.,thestationarysolutionofequation(3).Thissolutioniseffectivelythestationarydistribu- tionoftheGSRstatistic,(Rtr)t(cid:62)0,underthe“no-drift”hypothesis;seealso,e.g.,(Feinberg andShiryaev2006,Remark4.7,pp.466–467)and(Burnaevetal2009,Remark2,p.529). For a proof of existence of this limiting invariant distribution for any initial point Rr = 0 r(cid:62)0,see,e.g.,(PollakandSiegmund1985,Proposition3,p.271).Alltheseresultsarere- viewedatgreaterlengthinSection2whichisintendedtosurveytherelevantpriorliterature. Remarkably,theKolmogorovforwardequation(3)(possiblysubjecttodifferentconditions andconstraints)hasariseninotherdisciplinesaswell,notablyinphysics(inparticular,in quantummechanics),chemistry,andinmathematicalfinance.Itshouldthereforecomeas no surprise that in these fields the respective solution has been obtained using completely different techniques (e.g., Feynman path integrals; see Feynman 1948) and with no refer- encetotheGSRstatistic.InSection2webrieflyreviewtheseresultsaswell,astheymay castlightonapossibleinterpretationoftheGSRstatistic. The plan for the remainder of the paper is as follows. First, in Section 2 we give a brief account of the relevant prior literature. The main section of the paper is Section 3, which is where we first set up the corresponding Kolmogorov forward equation (3) more formally,andthensolveitexplicitly,i.e.,obtainanalyticallyaclosed-formformulaforthe P -transitionpdfoftheGSRstatistic.Next,inSection4theobtainedformulaisreconciled ∞ with several previously published parallel results. In Section 5 we exploit the found pdf p(x,t|r)numericallyandofferabriefnumericalstudyofthestatisticalbehaviorexhibited bytheGSRstatisticinthe“drift-free”regime.Tocarryoutthestudy,weimplementedthe obtainedpdfformulainMathematica1.Lastly,Section6ismeanttodrawalineunderthe entirepaper. 2 OverviewoftheRelevantPriorLiterature As can be gathered from the introduction, the centerpiece of this work is the solution of the Kolmogorov forward equation (3). This problem—in one or another form, shape, and 1 MathematicaisapopularsoftwarepackagedevelopedbyWolframResearch,Inc.asaprogrammingen- vironmentforscientificcomputing.SeeontheWebathttp://www.wolfram.com/mathematica/. 6 PolunchenkoandSokolov context—hasbeenaddressedintheliteraturebefore,andthissection’sobjectiveistoprovide abriefoverviewoftherelevantpriorresults. One of the first to arrive at a closely related problem was Wong (1964). The actual objective was to find the transition pdf of a stochastic process with a given stationary (as t → +∞) distribution. Specifically, the point of departure for Wong (1964) was the celebrated work of Kolmogoroff (1931), where it was noticed for the first time that it is possible to construct a diffusion process whose stationary distribution is a member of the Pearson (1895) distribution family2. However, since Kolmogoroff (1931) didn’t carry out theactualconstruction,Wong(1964)effectivelypickedupwhereKolmogoroff(1931)left off, and built on to the work of Karlin and McGregor (1960) to derive a number of pro- cesseswithaPearson-typestationarydistribution;foranoverviewofanalyticallytreatable Kolmogorov–Pearsondiffusionssee,e.g.,Avrametal(2012).OneofthePearsondistribu- tions considered by Wong (1964) was the inverted Gamma distribution, which is a gener- alizedextremevalueFréchet–Gumbel-typedistribution.Incidentally,thislatterdistribution wasdemonstratedbyShiryaev(1961,1963)tobethestationarydistributionoftheclassical SR statistic (special case of the GSR statistic with no headstart) set up for the Brownian motion change-point scenario (1) when the observed process, (Xt)t(cid:62)0, is still drift-free; seealso,e.g.,(FeinbergandShiryaev2006,Remark4.7,pp.466–467)and(Burnaevetal 2009,Remark2,p.529).Foraproofofexistenceofthisdistributionsee,e.g.,(Pollakand Siegmund1985,Proposition3,p.271)and(Peskir2006,Section2.2,p.3).Hence,effec- tively (Wong 1964, Section 2.F, p. 271) inadvertently discovered the SDE whose solution (inaspecialcase)istheGSRdiffusion3.One“problem”withWong’s(1964)work,how- ever, is that it is rather concise; in fact, it is so brief that it is, in essence, a “cook book” of“ready-made”formulaewithpracticallynoderivationthereofoffered.Thatsaid,Wong (1964)doesprovidethegeneralrecipetofindthetransitionpdfformulae,andtherecipeis tosolvetherespectiveKolmogorovforwardequationdirectlythroughtheFourierspectral methodwhichwementionedemergedfromtheWeyl–Titchmarshtheory;seeWeyl(1910) andTitchmarsh(1962).ItisthisapproachandmethodthatweintendtoemployinSection3 below, but with no key details omitted. Another “issue” with Wong’s (1964) work is that littleattentionispaidtotheamenabilityoftheobtainedtransitionpdfformulaetonumerical evaluation,anaspectimportantforapplications.Weprovideashortnumericalstudyofthe founddistributioncarriedoutwiththeaidofaMathematicascriptweprepared. Amorerecentandmorerelevantexampleofpriorworkonourproblemwouldbethe workofPeskir(2006).Justlikeusinthepresentpaper,Peskir(2006)alsofocusedon(effec- 2 Tobemoreprecise,Pearson(1895)onlyobtainedanddealtwithaspecialcaseofwhatlaterbecame knownasthePearson(first-orderseparabledifferential)equation,thebasisforthePearsondistributionfamily. 3 Wealsowarnofatypoonpage271of(Wong1964,Section2.F,p.271)inthesecondformulafromthe bottomofthepage(theformulaisunnumbered).Specifically,asgivenby(Wong1964,Section2.F,p.271), theformula(intheoriginalnotation)is (cid:18) α 1(cid:19) 2F0(−α−iµ,−α+iµ,−x)=xα+iµΨ − −iµ,1−2iµ, , 2 x √ wherei(cid:44) −1,2F0isthegeneralizedhypergeometricseries,andΨistheTricomiconfluenthypergeomet- ricfunction.ThetypoisthatthefirstargumentoftheΨ-functionintheright-handsideshouldbe−α−iµ, not−α/2−iµ.Thatis,thecorrectidentityis (cid:18) 1(cid:19) 2F0(−α−iµ,−α+iµ,−x)=xα+iµΨ −α−iµ,1−2iµ, , x ascanbeconfirmed,e.g.,with(AbramowitzandStegun1964,Section13.1,p.504). ExactDistributionoftheGeneralizedShiryaev–RobertsDiffusion 7 tively)theGSRstatistic,andconsideredaproblemevenmoregeneralthantheonetreated inthiswork.Specifically,Peskir(2006)wasafterthefundamentalsolutionofthefollowing Kolmogorovforwardequation: ∂ p˜(x,t|y)=− ∂ (cid:2)(1+ax)p˜(x,t|y)(cid:3)+b ∂2 (cid:2)x2p˜(x,t|y)(cid:3), t(cid:62)0, (5) ∂t ∂x ∂x2 wherea∈Randb(cid:62)0aregivenconstants.Thisequationismoregeneralthanequation(3) thatweareafter:thetwocoincidewhena=0andb=µ2/2.However,evenwhena(cid:54)=0, theforegoingequationisstillcloselyconnectedtotheGSRdiffusion:ifa = θµ,thenthe equationgovernsthetransitionpdfoftheGSRstatistic,(Rtr)t(cid:62)0,whenthedriftoftheob- servedprocess(Xt)t(cid:62)0 isθforallt (cid:62) 0,buttheGSRstatisticissetupto“anticipate”the drifttobeµ.Moreover,Peskir(2006)didnotrestrictattentiontothecasewhenx,y (cid:62)0(as wedointhiswork),butallowedx,y ∈R.OnecanthereforearguethatPeskir(2006)was, in effect, concerned with the transition pdf of the GSR statistic (with the headstart given and arbitrary from the entire real line) in a regime general enough to include as special cases(a)thepre-changeregime,(b)thepost-changeregime,and(c)thepost-changeregime withdrift-shift-misspecification,i.e.,whentheactualdrift-shiftmagnitudeisdifferentfrom theone“anticipated”bytheGSRstatistic;thislatterinterpretationcanbegivenbasedon, e.g., the work of Pollak and Siegmund (1985), where the performance of the original SR procedurewasstudiedintheminimaxBrownianmotioncontextwithapossibledrift-shift- misspecification. Peskir’s (2006) approach to solving equation (5) involved obtaining the Laplacetransform(theGreen’sfunction)oftheequation’sfundamentalsolution,andthen invertingittogetthesoughttransitionpdfitself.However,althoughPeskir(2006)didsuc- ceedinobtainingtheLaplacetransform(whichinthepre-changeregimetunedouttobethe LaplacetransformoftheHartman–Watsondistribution;seeHartmanandWatson1974),he wasabletoinvertitexplicitlyinthepre-changeregimeonlyandonlyforthezeroheadstart (i.e.,fortheclassicalSRstatistic).AsperthehopeexpressedbyPeskir(2006)forsomeone to“runthelastleg”andcompletehiswork,wegeneralizetheresultofPeskir(2006)forthe pre-changeregimetoanarbitrary(nonnegative)headstart. WenowrecalltheremarkwemadeintheintroductionthattheverysameKolmogorov forward equation (3) that we are after in this paper has also appeared in such disciplines asphysicsandmathematicalfinance.Specifically,inphysics,theequationisknownasthe Fokker–Planckequation(afterFokker1914andPlanck1917),andthereisvoluminouslit- eraturedevotedtothesolutionofthisequation(usingcompletelydifferentmethods).Tobe more specific, equation (3) arises in physics (as well as in chemistry), e.g., in connection withtheso-calledMorse(1929)potential.Thelatterisaharmonicoscillatormodelwidely usedtodescribethevariationofenergywithrespecttotheinternucleardistanceinadiatomic molecule.See,e.g.,Grosche(1988)andMorseandFeshbach(1953).Morespecifically,the respective equation that governs these energy variations is equivalent to our Kolmogorov forward equation (3) up to the initial and boundary conditions. It is known as one of the veryfewequationsinquantummechanicsthatpermitanalyticalsolution.Thestandardso- lutionstrategyusedinphysicstohandlethisequation(aswellasFokker–Plankequations ingeneral)istoexploittheFeynman(1948)pathintegralframework.Onahigherlevel,the physicalequivalentofourKolmogorovforwardequation(3)describesthepositionofapar- ticlemovingaroundinaninhomogeneousenvironment,drivenbyacombinationofrandom forces(e.g.,thermalnoise).Inthiscontext,equation(3)hasbeenstudied,e.g.,byMonthus andComtet(1994)andbyComtetandMonthus(1996). Another area to which the Kolmogorov forward equation (3) is no stranger is mathe- maticalfinance.Thefinancialinterpretationoftheequationisthatitdescribesthepriceof 8 PolunchenkoandSokolov acertainfinancialderivativeasafunctionoftimeandthederivativeparameters.Moreover, Markovprocesses(especiallyindimensionone)areoftenusedinfinancetomodelthetime- evolutionofthe“stateoftheeconomy”;thelatterisessentiallythefinancialequivalentofthe randomenvironmentinsidewhichthephysicalparticleismoving,ifoneistorelatetothe physicalanalogy.TheWeyl–TitchmarshspectraltheoryaswellastheFeynmanpathinte- gralapproachhavebothbeenusedinmathematicalfinancetopricefinancialinstrumentsand to describe the evolution of the underlying economic state. Excellent references assessing thestate-of-the-artoffinancialderivativepricingmethodswouldbeBorodinandSalminen (2002),Linetsky(2006,2007),andthereferencestherein.Bywayofexample,theequation consideredbyPeskir(2006),i.e.,equation(5),arisesinfinanceinthecontextofpricingthe so-calledAsianoptionsandMerton’s cashdividends;see,e.g.,Linetsky(2004a,b).Asan economicmodel,equationsakinto(3)orto(5)ariseinconnectiontotheso-calledHull– Whitevolatilitymodel(seeHullandWhite1987)whichpillarsonthesametypeofSDEs that governs the GSR statistic. For a thorough mathematical treatment of the Hull–White volatilitymodel,see,e.g.,Fatoneetal(2013).Moreover,diffusionprocessessimilartothe GSRstatistichavealsobeenencounteredasstochasticinterestratemodels,andhavebeen used to find the present value of annuities. See, e.g., Vanneste et al (1994); De Schepper et al (1994); De Schepper and Goovaerts (1999). In the latter papers, the authors devised theFeynmanpathintegralapproachtosolvetheequivalentofequation(3).Wenotethatin Feynmanpathintegralframeworkthesolutionneednotintegratetoone,asisrequiredinour casebythenormalizationconstraint(4).Asaresult,theobtainedsolutionsimplydiesoffto zeroastimegoeson,andnonontrivialstationarydistributionisexhibited,whichisinstark contrastwiththeaforementionedresultofShiryaev(1961,1963)concerningthe(nontrival) stationarydistributionoftheGSRstatistic. One more area where the Kolmogorov forward equation (3) has appeared is stochas- tic processes. In the latter field, one of the key methods to deal with the equation is the Feynman–Kac formula that emerged from the pioneering work of Feynman (1948) and that of Kac (1949, 1951, 1959). The Feynman–Kac formula establishes a bridge between parabolicpartialdifferentialequations(PDEs)andstochasticprocesses;wenotethatequa- tion(3)isaparabolicPDE.Anotherpossibleapproachtofindingthetransitionpdfofthe GSRdiffusion,(Rtr)t(cid:62)0,istoderiveitdirectlyfromthedefinitionoftheprocess,viz.from formula(2).Specifically,theideaistoexploitthelinearconnectionbetween(Rtr)t(cid:62)0 and thelikelihoodratioprocess{Λt}t(cid:62)0,andthencapitalizeonthefactthatΛt isageometric (exponential)Brownianmotion,awell-studiedtypeofdiffusions.Theaforementionedwork ofPeskir(2006)containsanumberofreferenceswherepreciselythisapproachwasdevised to derive the distribution of effectively the GSR diffusion (although for obvious reasons no connection to the GSR statistic was made). See, e.g., Yor (1992), Donati-Martin et al (2001), Dufresne (2001) and Schröder (2003). However, the obtained distribution is only forthecaseofnoheadstart,i.e.,assumingthestartingpointofthediffusioniszero. Weconcludethissectionwitharemarkthattheaboveisjustasmallsampleofexamples ofcontextswheretheKolmogorovforwardequationarises:the“footprint”oftheequation ismuchlargerandspansmanymoredisciplines.Itwouldbeadauntingtasktoreviewallof therelevantapplicationsoftheequationinaholisticmannerinasinglepaper. ExactDistributionoftheGeneralizedShiryaev–RobertsDiffusion 9 3 TheMainResult Thissectionisthecenterpieceofthiswork.Itisintendedtoprovideasolutiontothemain problemofthispaper:toobtainaclosed-formformulafortheP -transitionpdfoftheGSR ∞ diffusion,(Rtr)t(cid:62)0,definedby(2)above. Togetstarted,recallthatunderprobabilitymeasureP∞,theGSRdiffusion,(Rtr)t(cid:62)0, solvestheSDE:dRr = dt+µRrdB ,t (cid:62) 0,whereRt = r (cid:62) 0.Sincethedriftfunc- t t t r tion a(x,t) ≡ 1 and the diffusion coefficient b(x,t) = µx are both independent of time t (cid:62) 0,theGSRdiffusion,(Rtr)t(cid:62)0,istime-homogeneous.Consequently,accordingtothe fundamentalworkofKolmogoroff(1931),thesoughttransitionpdf,p(x,t|r),x,r,t (cid:62) 0, satisfiestheKolmogorovforwardequation: ∂ p(x,t|r)=− ∂ p(x,t|r)+ µ2 ∂2 (cid:2)x2p(x,t|r)(cid:3), x,r,t(cid:62)0, (6) ∂t ∂x 2 ∂x2 subjectto(a)the(natural)normalizationconstraint (cid:90) ∞ p(x,t|r)dx=1, r,t(cid:62)0, (7) 0 and(b)theinitialconditionlim p(x,t|r)=δ(x−r)validforallx.Sinceequation(6) t→0+ isaPDEoforderoneintimetandordertwoinspacex,andistoholdovertheentirespace (x,t)∈[0,∞)×[0,∞),thetemporalinitialconditionandthenormalizationconstraintare tobecomplementedbytwospacialboundaryconditions:oneatx=0(orasx→0+)and oneasx→+∞.Toobtaintheseconditionslet J(x,t|r)(cid:44)p(x,t|r)− µ2 ∂ (cid:2)x2p(x,t|r)(cid:3) 2 ∂x betheprobabilitycurrent,sothatequation(6)canberewrittenmorecompactlyas ∂ ∂ p(x,t|r)+ J(x,t|r)=0, (8) ∂t ∂x whichiseffectivelya(one-dimensional)continuityequation:itconstitutesthelawofcon- servationofprobability,ananalogueofthewell-knownlawofconservationofenergyfrom physics.Forfurthernotationalbrevity,weshallomittheheadstartrinp(x,t|r)andsimply writep(x,t)throughoutthissection,unlessitisnecessarytoemphasizethedependenceon r. Toobtainthefirstboundarycondition,observethatsinceRr =r(cid:62)0byassumption,it 0 canbededucedfromdefinition(2)thatRr (cid:62)0almostsurelyunderanyprobabilitymeasure t andforanyt (cid:62) 0.IfRr isinterpretedasthepositionofahypothetical“particle”attime t instancet,thenthenonnegativityofRr forallt (cid:62) 0istosaythattheparticleistonever t leave the nonnegative half-plane. Put another way, no particle flow through the x-axis is permitted.Asaresult,theprobabilitycurrentJ(x,t)throughthex-axismustbezeroatall times,i.e.,J(0,t)=0forallt(cid:62)0,orexplicitly p(x,t)− µ2 ∂ (cid:2)x2p(x,t)(cid:3)=0, x→0+, (9) 2 ∂x whichisthefirstboundarycondition. The second boundary condition can be obtained by integrating both sides of (8) with respecttoxover[0,+∞)andthenusingthenormalizationconstraint(7).Specifically,this 10 PolunchenkoandSokolov yieldsthatJ(+∞,t)−J(0,t)=0forallt(cid:62)0,.Therefore,theprobabilitycurrentthrough x = +∞mustmatchtheprobabilitycurrentthroughx = 0atalltimes.Sinceinthelatter casetheprobabilitycurrentiszero,thesecondboundaryconditionis: p(x,t)− µ2 ∂ (cid:2)x2p(x,t)(cid:3)=0, x→+∞. (10) 2 ∂x Tosolveequation(6),letusfirstconsidertheequation’sstationarysolution.Specifically, thestationarityhereisinthetemporalsense,i.e.,p(x,t)isindependentoft,whichisthe caseinthelimitast→+∞.Letρ(x)(cid:44)lim p(x,t)denotethecorrespondinglimit. t→+∞ Foraproofthatρ(x)exists,see,e.g.,(PollakandSiegmund1985,Proposition3,p.271). Sinceinthestationaryregime∂p(x,t)/∂t ≡ ∂ρ(x)/∂t ≡ 0forallx,equation(8)simpli- fiesto∂J(x,t)/∂x = 0.Thelatterequation,inturn,istosaythattheprobabilitycurrent, J(x,t),asafunctionoftime,isconstant,andinviewofthetwoboundaryconditionsestab- lishedabove,thevalueofthatconstantiszero.Hence,weobtain ρ(x)− µ2 d (cid:2)x2ρ(x)(cid:3)=0. 2 dx Thesolutiontotheforegoingequationiswell-knownandisgivenbytheFréchet-type distributiondensity ρ(x)=e−µ22xµ22x2 1l{x(cid:62)0}. (11) ThisresultwasfirstobtainedbyShiryaev(1961,1963)asthestationarydistributionof theoriginalSRstatistic(withzeroheadstart).Seealso,e.g.,FeinbergandShiryaev(2006) and(Burnaevetal2009,Remark2,p.529).Asastationarydistribution,notonlyisρ(x) independentoft,butitisalsoindependentoftheheadstartr.Therefore,ρ(x)givenby(11) is the stationary distribution of the GSR statistic, (Rtr)t(cid:62)0, as well. It is also noteworthy thatthedistribution(11)isaspecialcaseoftheextremevalueinverseGammadistribution, andinacompletelydifferentcontext(viz.finance)andusingdifferenttechniqueshasbeen discovered,e.g.,by(Milevsky1997,Theorem1,p.224).Itisalsoaspecialcaseof(Comtet etal1998,Formula(74),p.264). Letusnowattackequation(6)forarbitrary(finite)time0 (cid:54) t < ∞.Tothatend,let usfirstheuristicallyoutlineourstrategy.Themainideaistosupposethatp(x,t)isofthe formp(x,t)=ρ(x)ψ(x)τ(t)whereρ(x)isthestationarydensitygivenby(11),andψ(x) andτ(t)aretobefound.Inotherwords,theideaistoassumethatthespacialvariable,x, and the temporal variable, t, can be separated. If that were the case, then the substitution p(x,t)=ρ(x)ψ(x)τ(t)wouldbringequation(6)intothefollowingform: τ(cid:48)(t) = 1 (cid:18)− d (cid:2)ρ(x)ψ(x)(cid:3)+ µ2 d2 (cid:2)x2ρ(x)ψ(x)(cid:3)(cid:19), x,t(cid:62)0. τ(t) ρ(x)ψ(x) dx 2 dx2 Sincethexandtvariablesarenowondifferentsidesoftheequation,inorderforthe twosidestobeequaltooneanotherirrespectiveofxandt,thetwosidesarebothtobeequal tothesameconstant,sayλ.Therefore,thesubstitutionp(x,t)=ρ(x)ψ(x)τ(t)effectively splitstheoriginalPDE(6)intotwoordinarydifferentialequations(ODEs): τ(cid:48)(t) =λ and 1 (cid:18)− d (cid:2)ρ(x)ψ(x)(cid:3)+ µ2 d2 (cid:2)x2ρ(x)ψ(x)(cid:3)(cid:19)=λ, (12) τ(t) ρ(x)ψ(x) dx 2 dx2