An Economic Bubble Model and Its First Passage Time
AAn Economic Bubble Model and Its First Passage Time
Angelos Dassios ∗ and Luting Li † Department of Statistics, London School of EconomicsMarch 23, 2018
Abstract
We introduce a new diffusion process { X t } t ≥ to describe asset prices within an economic bubblecycle. The main feature of the process, which differs from existing models, is the drift term wherea mean-reversion is taken based on an exponential decay of the scaled price. Our study shows thescaling factor on { X t } t ≥ is crucial for modelling economic bubbles as it mitigates the dependencestructure between the price and parameters in the model. We prove both the process and its firstpassage time are well-defined. An efficient calibration scheme, together with the probability densityfunction for the process are given. Moreover, by employing the perturbation technique, we deduce theclosed-form density for the downward first passage time, which therefore can be used in estimatingthe burst time of an economic bubble. The object of this study is to understand the asset pricedynamics when a financial bubble is believed to form, and correspondingly provide estimates to thebubble’s crash time. Calibration examples on the US dot-com bubble and the 2007 Chinese stockmarket crash verify the effectiveness of the model itself. The example on BitCoin prediction confirmsthat we can provide meaningful estimate on the downward probability for asset prices. Keywords : Economic Bubbles, Diffusion Process, First Passage Time, Perturbation, Cryptocurrency An economic bubble usually refers to economic phenomenons that asset prices extremely deviate fromtheir fundamental values [38]. One of the most famous bubbles in history, known as the Dutch TulipBubble [12, 18], could be traced back to the 1630s. According to P.M. Garber [18], from November 1636to February 1637, the prices of tulip bulbs had increased about 20 times. At the peak of the bubble, byselling a few bulbs people could even buy a luxury house in Amsterdam. However, only three monthslater, the bulbs became worthless. The rapid increases and sudden drops in asset prices are a commonfeature reflected by a bubble cycle. More modern examples can be found in [44, 25, 22].The burst of an economic bubble sometimes follows with financial crisis, or even economic depression.In modern history, the most devastating crisis would be the 2007-2009 Financial Crisis [36], where people ∗ [email protected] † [email protected] a r X i v : . [ q -f i n . M F ] M a r elieve the crash of the US real estate market is one of the causing. And the crash itself, is usuallyreferred to as the burst of the US Housing Bubble [22]. Although it is believed that a bubble cannot bepredicted before it is formed, by knowing the burst time in advance, governments and market participantscan manage the potential risk accordingly. Therefore, an effective estimate before the crash will help inpreventing systematic risk. The object of this study is to understand the asset price dynamics within aneconomic bubble cycle and provide estimates to the probability distribution of the collapse time.The financial bubbles have been studied extensively in econometrics and statistics. As a non-conclusivereview, we refer to [35] and the literatures it mentioned for the econometric approach; agent-based modelsin statistics and a summary of literatures can be found in [16]. In financial mathematics, local martingalemodels have been considered in option pricing problems. A. Cox and D.G. Hobson [10] included a widebranch of stochastic diffusions in their work. S. Heston et al. [20] enriched the discussions by introducingCIR process and Heston stochastic volatility model. In terms of the burst time prediction, C. Brooks andA. Katsaris [7, 6] forecasted the collapse of speculative bubbles in S&P 500 index using a three-regimemodel. To the best of our knowledge, there is limited research in modelling economic bubble dynamicsvia a pure time-homogeneous diffusion process. The research on finding the explicit probability densityof bubble crash time is even less. One paper related to our work is contributed by A. Kiselev and L.Ryzhik [27], where a mean-reversion process with an exogenous functional drift has been considered.In this paper we introduce a new time-homogeneous diffusion model. Our motivation is to providean alternative approach, where with the mathematical form to be as simple and tractable as possible, toprobabilistically describe asset price dynamics within an economic bubble cycle. The new model is closelylinked to the Shiryaev process [40, 39] derived by A.N. Shiryaev in the context of sequential analysis. Ourmodel involves three independent parameters, where two of them provide mean-reversion effects as in theOrnstein-Uhlenbeck (OU) process. As a crucial variable to our model, the third parameter controls thespeed of exponential decay in the drift term. Consequently, the dependence structure among the return,asset price, equilibrium level and the mean-reversion rate has been mitigated. Without introducing extrafunctionals, our model provides sufficient degree of freedom for calibrations, while on the other hand,avoids over-fitting. Due to the simple structure of the model, we are able to show the closed-form densityfunction of the bubble crash time.The main contribution in the present paper is that we have provided a self-contained material inmodelling bubble dynamic and predicting its burst time. On the theoretical side, we have proved the newmodel is a well-defined diffusion process and its first passage time (FPT) exists. To be more specific, theprocess is a semimartingale with a strong and unique solution. As a recurrent strong Markov process, themodel embeds an a.s. finite FPT; and its stationary distribution has been found with a neat functionalform. On the practical side, a calibration algorithm based on economic features has been considered.We have given explicit solution to the distribution of the process at fixed time. Moreover, the Laplacetransform (LT) of the FPT has been found, and based on the perturbation technique we have solved theclosed-form density for the downward FPT. In the end, the effectiveness of the model and its FPT density(FPTD) has been verified by three numerical examples.The rest of the paper is organised as follows: Section 2 introduces the SDE of our new model and themotivation behind it; Section 3 discusses the theoretical results from the new process itself; the closed-2orm solution of the FPTD is given by Section 4; in Section 5 we demonstrate the calibration algorithmand illustrate the model application via three examples, among which a prediction on the BitCoin collapsetime has been given; Section 6 concludes. Consider a filtered probability space { Ω , F , P } , where F = {F t } t ≥ is a natural filtration generated by astandard Brownian motion { W t } t ≥ . We introduce the following three-parameter SDE dX t = (cid:15) (cid:0) e − αX t − c (cid:1) dt + dW t , X = x ∈ R . (1)The parameters (cid:15), α are restricted on the positive real line and 0 ≤ c ≤ α controls the speed, curvature, and higher order information inthe drift term. Figure 1 illustrates the functionals of e − αX t to different choices of α . We can see, as afunction of X t , small α produces mildly linear decays in the drift. This extends the range of the processwhere positive return is maintained. On the other hand, large α generates evident exponential decays.In this case the drift sign is sensitive to the values of X t , and the range of positive drift is compressed.Figure 1: Functoin plots of e − αX t with α =0 . , . , ,
2. Green zone: positive drift; red zone:negative drift. Figure 2: Sample path for X t in 4 years time. Pa-rameters are chosen as α = 1, (cid:15) = 0 . c = 0 . X = 0 and dt = .To illustrate the new process in a more intuitive way, we plot the simulated path of X t with α = 1(green curve in Figure 1) in 4 years time. The other parameters are chosen as (cid:15) = 0 . c = 0 . X = 0.Three thresholds in colors of (from below to above) green, black and red indicate different regimes forthe process: I) when X t is negative or near 0 (green line), according to SDE (1) the process embeds astrong positive trend; II) black line plots the equilibrium level where e − αX t = c and in a long-run X t oscillates around this position; III) red line shows the level of X t where e − αX t = 0 . c and the process3s forced to drop back due to the strong negative trend. In addition, the sample path shows X t spendsmuch less time in visiting the equilibrium level from the initial point, than that it drops back from thesymmetrical high position. Consider the green curve in Figure 1, this asymmetric feature is a naturalreflection to the exponential transforms from the level X t to the instantaneous return. As a result, X t ingeneral should have a rapid increase when it is below the equilibrium level, but, even though X t exceedsthe equilibrium level to higher positions, it is not necessary that X t will drop down immediately. Thisbehaviour essentially differentiates our new model and the OU type mean-reversion processes.The model feature coincides with observations from economic bubbles. Refer to the theory by H.P.Minsky and H. Kaufman [30]. A bubble cycle is formed by five steps: Displacement, Boom, Euphoria,Profit Taking, Panic . In the first step the asset price remains at a lower level and the process usuallyhas an ‘initiative’ to increase. This corresponds to the regime I in our model. During the booming stage,the price becomes sensitive to positive market news and increases rapidly. Although sometimes due todivergence in market anticipations that the price may drop down, after oscillations the asset price willkeep increasing. This is described by regime II. In regime III, the peak is shown and large negativedrift is accumulated. This describes the euphoria stage where the asset price hits historically high levels;however, due to market capital limits or aversions of risk, the market expectations become negative. Inthe profit taking stage, asset price becomes sensitive to negative market news and the process shifts fromregime III to regime II. In the end the process drops back to the mean-reversion level, or even continueto drop to regime I. This describes the last step of the bubble.From the calibration point of view, α mitigates the dependence structure of the instantaneous returnto the price, equilibrium level and mean-reversion rate. Consider the drift function where α is suppressed, µ ( X t ) := (cid:15) ( e − X t − c ) . In this model once c is determined, the equilibrium level X t = − ln( c ) becomes a fixed number. If a largerate of c is calibrated, then we simultaneously have a small equilibrium level − ln( c ). Therefore when X t is small, where e − X t − c is close to 0, in order to fit a large instantaneous return the mean-reversion rate (cid:15) should be adjusted highly as well. But we know usually a bubble spends years to finish its whole cycle;so large reversion rate is not desired for a bubble model.The analysis shows the α -suppressed model is not capable for calibrating a bubble dynamic. As acomplement, extra functional term is required (cf. [27]). However, without introducing extra functionsour three-parameter model extends the freedom in model calibration. Combining previous discussions wesee SDE (1) is a good candidate for describing economic bubbles. Otherwise the process will take a long period to visit regime III), where e − X t ≈ − c becomes dominating. Theoretical Soundness
Proposition 3.1
There exists a unique and strong solution { X t } t ≥ to SDE (1) and which has thefollowing explicit form X t = x + W t − c(cid:15)t + 12 α ln (cid:18) (cid:15)αe − αx (cid:90) t e − α ( W s − c(cid:15)s ) ds (cid:19) ; moreover { X t } t ≥ is a strong Markov process. Proof:
We consider an exponential transform on X t such that Y t = e αX t with Y = e αx . By applying Ito lemma we show the coefficients of Y t satisfy global Lipschitz continuity and linear growthconditions: dY t = 2 α [ (cid:15) + ( α − c(cid:15) ) Y t ] dt + 2 αY t dW t . (2)According to Theorem 2.9 in [26] we conclude there exists a unique and strong solution { Y t } t ≥ to SDE(2). Therefore { X t } t ≥ is the unique and strong solution to SDE (1). On the other hand, refer to [43,Section 4.4] Y t has the following explicit form: Y t = e α ( W t − c(cid:15)t ) (cid:20) Y + 2 α(cid:15) (cid:90) t e − α ( W s − c(cid:15)s ) ds (cid:21) . Then by substituting Y t = e αX t into the equation above we solve X t .Now we consider the strong Markov property. Note that the coefficients in SDE (2) are continuousso are bounded on compact subsets of R . Combining the well-posed proof in above, and referring to [26,Theorem 4.20] we show { Y t } t ≥ is a strong Markov process. Therefore the strong Markov property holdsfor { X t } t ≥ . (cid:3) Remark 3.2
The proof depicts { X t } t ≥ from another aspect. SDE (2) shows { Y t } t ≥ is a geometricBrownian motion with a mean-reversion drift. Referring to [40, Equation (9)] this is indeed a Shiryaevprocess. Therefore { X t } t ≥ is the logarithm of the Shiryaev process. From the explicit solution in Proposition 3.1 we see { X t } t ≥ is a semimartingale, where the boundedvariation (BV) part consists of a strictly decreasing function and a strictly increasing function. Dependingon the Brownian motion path in the exponential integral, for different t > c = 0 it is clear that only the increasing function isretained. This indicates under special circumstance { X t } t ≥ could be a submartingale. Corollary 3.3 If c = 0 then { X t } t ≥ is a strict submartingale. roof: When c is suppressed the solution of { X t } t ≥ becomes X t = x + W t + 12 α ln (cid:18) (cid:15)αe − αx (cid:90) t e − αW s ds (cid:19) . (3)The adeptness is clear from definition. We consider the L -integrability of X t . Note that by applying theJensen’s inequality for concave function ln( · ), we have E (cid:20) ln (cid:18) (cid:15)αe − αx (cid:90) t e − αW s ds (cid:19)(cid:21) ≤ ln (cid:18) E (cid:20) (cid:15)αe − αx (cid:90) t e − αW s ds (cid:21)(cid:19) . By changing integral and expectation, E (cid:20)(cid:90) t e − αW s ds (cid:21) = 12 α (cid:16) e α t − (cid:17) . (4)On the other hand, when α, (cid:15) ∈ (0 , + ∞ )1 + 2 (cid:15)αe − αx (cid:90) t e − αW s ds ≥ , ∀ t ≥ . (5)Therefore (cid:12)(cid:12)(cid:12)(cid:12) ln (cid:18) (cid:15)αe − αx (cid:90) t e − αW s ds (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) = ln (cid:18) (cid:15)αe − αx (cid:90) t e − αW s ds (cid:19) . (6)Combining (4) and (6) we have12 α E (cid:20)(cid:12)(cid:12)(cid:12)(cid:12) ln (cid:18) (cid:15)αe − αx (cid:90) t e − αW s ds (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:21) ≤ α ln (cid:16) (cid:15)α e − αx (cid:16) e α t − (cid:17)(cid:17) . (7)In the end, note that E [ | W t | ] = (cid:114) tπ . (8)So applying the triangle inequality and combining (7) and (8) we show the L -integrability of X t by E [ | X t | ] ≤ x + (cid:114) tπ + 12 α ln (cid:16) (cid:15)α e − αx (cid:16) e α t − (cid:17)(cid:17) < + ∞ . Finally, the non-decreasing conditional expectation of E (cid:20) X t (cid:12)(cid:12)(cid:12)(cid:12) F s (cid:21) , for 0 ≤ s < t < + ∞ , is given by againusing (5), that ln (cid:18) (cid:15)αe − αx (cid:90) t e − αW s ds (cid:19) > . This concludes our proof. (cid:3) { X t } t ≥ We consider finding the distribution of { X t } t ≥ . By Proposition 3.1 we see the solution of X t involvesBrownian motion and its exponential integral. Similar problem has been answered by A. Dassios andJ. Nagaradjasarma [13] for the square-root process. G. Peskir [33] deduced the fixed time distributionfor the Shiryaev process in a special case. But for the general case only the Laplace transform has beengiven. Here we refer to the results about Brownian motion and its exponential integral in H. Matsumotoand M. Yor [29, 45], and have the following proposition.6 roposition 3.4 For fixed t > and u ∈ R , the probability density of X t is given by P ( X t ∈ du ) = αdu · (cid:20)(cid:90) ∞ ζ ( u ; c(cid:15)α , y ) exp (cid:18) − c (cid:15) t + 1 /y + ζ ( u ; 2 , y )2 (cid:19) θ (cid:0) ζ ( u ; 1 , y ) , α t (cid:1) dy (cid:21) , where θ ( r, s ) = r √ π s e π s (cid:90) ∞ e − v s − r cosh( v ) sinh( v ) sin (cid:16) πvs (cid:17) dv and ζ ( u ; µ, y ) = (cid:0) (cid:15)α e − αx y (cid:1) µ e − µ · α ( u − x ) y . Proof:
Let s = α t , then for another standard Brownian motion B s , with probability 1 the followingequation holds true B s + c(cid:15)α s = − α ( W t − c(cid:15)t ) . Denote by µ := c(cid:15)α , B ( µ ) s := B s + µs, A ( µ ) s := (cid:90) s e B ( µ ) v dv. (9)Then referring to [29, 45] we have P (cid:16) A ( µ ) s ∈ dy, B ( µ ) s ∈ dz (cid:17) = 1 y exp (cid:18) µz − µ s − e z y (cid:19) · θ (cid:18) e z y , s (cid:19) dydz, (10)where θ ( r, ξ ) = r (cid:112) π ξ e π ξ (cid:90) ∞ e − v ξ − r cosh( v ) sinh( v ) sin (cid:18) πvξ (cid:19) dv. On the other hand, re-express X t using B ( µ ) s and A ( µ ) s . Note that A ( µ ) s = (cid:90) α t exp (cid:16) − α (cid:16) W vα − c(cid:15) vα (cid:17)(cid:17) dv. By changing variable with w = vα we have A ( µ ) s = α (cid:90) t e − α ( W w − c(cid:15)w ) dw. (11)Rewrite { X t } t ≥ in Proposition 3.1 using (9) and (11), we get X t P = x − B ( µ ) s α + 12 α ln (cid:16) (cid:15)α e − αx A ( µ ) s (cid:17) . (12)We now consider the density function for X t . Note that for fixed u ∈ R and X t ≤ u , (12) implies α ( x − u ) + 12 ln (cid:16) (cid:15)α e − αx A ( µ ) s (cid:17) ≤ B ( µ ) s . Denote by g (cid:16) u, A ( µ ) s (cid:17) := α ( x − u ) + 12 ln (cid:16) (cid:15)α e − αx A ( µ ) s (cid:17) . Considering (10) we have P ( X t ≤ u ) = (cid:90) y ∈ (0 , ∞ ) (cid:90) z ≥ g ( u,y ) P (cid:16) A ( µ ) s ∈ dy, B ( µ ) s ∈ dz (cid:17) . (13)7aking derivatives on u we further get P ( X t ∈ du ) = − (cid:90) y ∈ (0 , ∞ ) y exp (cid:18) µg ( u, y ) − µ s − e g ( u,y ) y (cid:19) · θ (cid:18) e g ( u,y ) y , s (cid:19) dy · g (cid:48) u ( u, y ) du (14)= α (cid:90) y ∈ (0 , ∞ ) y exp (cid:18) µg ( u, y ) − µ s − e g ( u,y ) y (cid:19) · θ (cid:18) e g ( u,y ) y , s (cid:19) dydu. (15)Introduce the function ζ ( u ; µ, y ) that ζ ( u ; µ, y ) := e µg ( u,y ) y = (cid:0) (cid:15)α e − αx y (cid:1) µ e − µ · α ( u − x ) y . Then rearranging (15) we have1 y exp (cid:18) µg ( u, y ) − µ s − e g ( u,y ) y (cid:19) · θ (cid:18) e g ( u,y ) y , s (cid:19) = ζ ( u ; µ, y ) exp (cid:18) − µ s − /y + ζ ( u ; 2 , y )2 (cid:19) · θ ( ζ ( u ; 1 , y ) , s ) . The proof is concluded by substituting µ = c(cid:15)α and s = α t into the equation in above. (cid:3) Remark 3.5
The function θ ( r, s ) is closely related to the study in Hartman-Watson distributions [19].As noticed by H. Matsumoto, M. Yor [29], and other researchers [3, 23], θ ( r, s ) is highly oscillating,especially for small s . Therefore it is not easy to compute the accurate values of the density. In practice it is more meaningful to provide the probability distribution function rather than thedensity function. This requires an extra integral on P ( X t ∈ du ). Considering the integral involved in θ ( r, s ), and the integral taking on θ ( r, s ), in total we need to compute three integrals for the distributionfunction P ( X t ≤ du ). A direct finite-difference scheme would therefore generate computational efficiencyissue. Instead, we consider computing the probability via Monte Carlo simulation.Referring to Equation (13) and the density function in Proposition 3.4, we have two choices in de-veloping the simulation algorithm. Based on (13) we could follow the acceptance-rejection approach byconsidering the relative positions between z and g ( u, y ). However, in the present paper we will concentrateon the direct sampling scheme by employing the explicit density function. Proposition 3.6
For fixed t > and u ∈ R , define m ( z, y ) = αζ ( z ; c(cid:15)α , y ) exp (cid:18) − c (cid:15) t + 1 /y + ζ ( z ; 2 , y )2 (cid:19) · ˆ θ (cid:0) ζ ( z ; 1 , y ) , α t (cid:1) , where ˆ θ ( r, s ) := r s e π s E (cid:20) V e − r cosh( V ) sinh( V ) sinc (cid:18) Vs (cid:19)(cid:21) with sinc( w ) := sin( wπ ) wπ and V ∼ N (0 , √ s ) . Then for two i.i.d. uniformly distributed random variables U and Y , the probability distribution of X t is given by P ( X t ≤ u ) = E (cid:34) m (cid:0) − U + u + 1 , Y − (cid:1) U Y (cid:35) . roof: First we show the identity between ˆ θ ( r, s ) and θ ( r, s ). Recall Proposition 3.4 that θ ( r, s ) = r √ π s e π s (cid:90) ∞ e − v s − r cosh( v ) sinh( v ) sin (cid:16) πvs (cid:17) dv. Rewriting the function we get θ ( r, s ) = rs e π s (cid:90) ∞ ve − r cosh( v ) sinh( v ) sinc (cid:16) vs (cid:17) · √ πs e − v s dv = r s e π s E (cid:20) V e − r cosh( V ) sinh( V ) sinc (cid:18) Vs (cid:19)(cid:21) =:ˆ θ ( r, s ) . Note that the second equation holds true is due to the fact( v sinh( v )) · e − r cosh( v ) · sinc (cid:16) vs (cid:17) is an even function.Now let m ( z, y ) to be defined as in Proposition 3.6. Based on the identity between θ ( r, s ) and ˆ θ ( r, s ),and referring to Proposition 3.4, we can write the probability distribution of X t as P ( X t ≤ u ) = (cid:90) u −∞ (cid:90) ∞ m ( z, y ) dydz. (16)Change variables that U = − z − u − , Y = 11 + y . Then re-expressing z, y by U, Y in (16) we have P ( X t ≤ u ) = − (cid:90) (cid:90) m (cid:18) − U + u + 1 , Y − (cid:19) dU dYU Y . By noticing the fact that uniform distribution has constant probability density dU = dY = 1 we concludethe proof. (cid:3) Remark 3.7
The main consideration of involving sinc( · ) is to reduce the oscillation effects from func-tion sin( · ) . From the numerical calculation point of view, the function, though cannot totally solve theoscillating issue, could mitigate the chaos to some extent. Proposition 3.8
For t ↑ + ∞ , the stationary distribution of X ∞ := lim t ↑ + ∞ X t is given by p ( x ) = (cid:20)(cid:90) y ∈ R w ( x ) w ( y ) dy (cid:21) − , where w ( x ) = exp (cid:110) (cid:15)α e − αx + 2 (cid:15)cx (cid:111) . roof: Consider the Fokker-Planck equation at t = + ∞ :12 p (cid:48)(cid:48) ( x ) − (cid:15) ( e − αx − c ) p (cid:48) ( x ) + 2 α(cid:15)e − αx p ( x ) = 0 . Define w ( x ) as in the proposition. Solving the ODE without boundary conditions we get p ( x ) = C (cid:82) x −∞ w ( y ) dy + C w ( x ) . Note that w ( y ) ↑ + ∞ and is dominated by a double exponential function when y ↓ −∞ . As w ( y ) ≥ R , so for x > −∞ the integral does not exist: (cid:90) x −∞ w ( y ) dy = + ∞ . In order to get a valid density function we therefore set C = 0. Determining C by the full-integrabilitycondition we conclude the proof. (cid:3) Remark 3.9
The stationary distribution is right-skewed due to the fact that the double exponential func-tion diverges much faster than the exponential function. In fixed income modelling, the double exponentialfunction is also used in term structure calibrations (cf. Nelson-Siegel model [21]).
In the later section we will deduce the probability density function of the FPT for { X t } t ≥ . Beforeconducting the calculations we show the existence of the FPT to any constant level a ∈ R . Proposition 3.10 { X t } t ≥ is a recurrent process on R . Proof:
Consider the substitution that Z t := e − αX t , Z = e − αx . (17)Applying the Ito lemma we get dZ t = 2 αZ t ( − (cid:15)Z t + c(cid:15) + α ) dt − αZ t dW t . (18)Based on Proposition 3.1 we know { Z t } t ≥ is a diffusion process with unique and strong solution; moreoverthe strong Markov property holds as well .Our construction indicates that { Z t } t ≥ only takes value on the positive half-plane. Additionally bychecking with (18) we see Z t = 0 is an absorbing bound. So consider I = (0 , + ∞ )as the domain of { Z t } t ≥ . We show { Z t } t ≥ is recurrent on I using the scale function as discussed in[26]. For any fixed parameter A ∈ I , we define s ( z ) := (cid:90) zA exp (cid:40) − (cid:90) ξA ( − (cid:15)ζ + c(cid:15) + α ) αζ dζ (cid:41) dξ, z ∈ I. (19) In fact one can even find the explicit solution of Z t by referring to the stochastic Verhulst equation in [43]. z ∈ I the nondegeneracy condition4 α z > z ∈ I , consider δ > z − δ ∈ I . Then the localintegrability condition is satisfied by showing (cid:90) z + δz − δ αζ | ( − (cid:15)ζ + c(cid:15) + α ) | α ζ dζ ≤ δ α ( z − δ ) + (cid:15)α δ + c(cid:15) + α α ln (cid:18) z + δz − δ (cid:19) < + ∞ . Therefore the scale function (19) is well defined.We now calculate the limit value of the scale function at boundaries of I . Rewrite the s ( z ) in (19) as s ( z ) = A c(cid:15)α e − A(cid:15)α (cid:90) zA e (cid:15)α ξ ξ c(cid:15)α dξ. (20)Let l + := 0 + and r − := ∞ − . Substituting the boundary values into (20) we get s ( l + ) = − A c(cid:15)α e − A(cid:15)α (cid:90) A e (cid:15)α ξ ξ c(cid:15)α dξ = −∞ , and s ( r − ) = A c(cid:15)α e − A(cid:15)α (cid:90) ∞ − A e (cid:15)α ξ ξ c(cid:15)α dξ = + ∞ . Therefore according to Proposition 5.22 in [26] we conclude that { Z t } t ≥ is recurrent on I . The recurrenceof { X t } t ≥ on R follows by (17). (cid:3) Corollary 3.11
For any a, X = x ∈ R , the first hitting time of { X t } t ≥ from x to a exists. Proof:
This directly follows from Proposition 3.10. (cid:3)
Remark 3.12
Note that in Corollary 3.11 there is no restriction on the direction of FPT. More specifi-cally, define τ ax ↑ := inf { t ≥ X t = a | x < a } , and τ ax ↓ := inf { t ≥ X t = a | x > a } to be the FPTs for from below and from above respectively. Then P (cid:0) τ ax ↑ < + ∞ (cid:12)(cid:12) X = x ) = 1 and P (cid:0) τ ax ↓ < + ∞| X = x (cid:1) = 1 . { X t } t ≥ { X t } t ≥ is a well-defined diffusion process. Therefore the corresponding in-finitesimal generator exists. Denote by C the collection of twice differentiable and continuous functionsdefined on R . For f ∈ C , the infinitesimal generator of { X t } t ≥ is given by A f ( x ) = (cid:15) ( e − αx − c ) f (cid:48) ( x ) + 12 f (cid:48)(cid:48) ( x ) . (21)11 .1 Dirichlet Problem Note by our settings the filtration F is continuous on both sides. So it is equivalent to consider the FPTto either an open or a closed set. W.l.o.g. , for a ∈ R we define D u := { x ∈ R : x > a } , D l := { x ∈ R : x < a } to be the domains of upper and lower regions to a . For notational convenience we denote by D to referto either D u or D l .Let ∂ D to be the set of boundaries of D . Then ∂ D u := { a, + ∞} , ∂ D l := { a, −∞} . The FPT of { X t } t ≥ with X = x ∈ D can be defined correspondingly as τ ax := inf { t ≥ X t ∈ ∂ D} . For a short notation we suppress x, a and write τ := τ ax . The existence of τ is given by Corollary 3.11. In this section we follow G. Peskir and A.N. Shiryaev[34] to deduce the Dirichlet type boundary value problem for the Laplace transform of τ . Consider anarbitrary sequence of well-defined stopping times { σ n } ≤ n ≤ + ∞ . Define σ := lim n ↑ + ∞ σ n . Note that { X t } t ≥ is a continuous process. Therefore { X t } t ≥ is continuous over all stopping times, i.e.lim n ↑ + ∞ X σ n = X σ . Moreover, by Proposition 3.1 { X t } t ≥ is a strong Markov process. For fixed β ≥
0, define f ( x ) = E x (cid:2) e − βτ (cid:3) , x ∈ D , (22)where E x [ · ] := E [ ·| X = x ] = E [ ·|F ]. Then refer to [34] f ( x, β ) is the unique solution to the followingODE A f ( x ) = βf ( x ) , x ∈ D (23)with Dirichlet type boundary conditions f ( ∂ D ) = (1 , T . (24) Remark 4.1
The result in above follows from the killed version of potential theory. Note that the strongMarkov property and continuity over all stopping times are crucial for representing the unique solutionfrom the Dirichlet problem by (22) . Remark 4.2
The boundary conditions (24) imply the boundness in f , though the set C does not requirethe boundness explicitly. Therefore by constructing a proper martingale we can use Feynman-Kac theoremto deduce similar conclusion. However, in order to employ the optimal sampling theorem the boundness-related properties in τ should be further demonstrated. .2 Direct Solution to the Dirichlet Problem Let m := √ c (cid:15) +2 β − (cid:15)c α ,n := √ c (cid:15) +2 β + αα ,ψ := (cid:15)e − αx α ,λ := x ( (cid:15)c − (cid:112) c (cid:15) + 2 β ) . (25)Refer to [2]. The solution of (23), by substituting (21) into, is given by f ( x ) = C e λ M ( m, n, ψ ) + C e λ U ( m, n, ψ ) , (26)where M ( m, n, ψ ) and U ( m, n, ψ ) are solutions to the Kummer’s equation [11] ψu (cid:48)(cid:48) ( ψ ) + ( n − ψ ) u (cid:48) ( ψ ) = au ( ψ ) . Now determine the constants C and C . Consider the hitting from below case, i.e. the boundary istaken on ∂ D l . Substituting x = −∞ we see ψ = + ∞ and λ = + ∞ as (cid:15)c − (cid:112) c (cid:15) + 2 β <
0. Refer to [28], the asymptotic of U ( m, n, ψ ) for large ψ is given by U ( m, n, ψ ) ∼ ψ − m , ψ ↑ + ∞ . Although by m > U ( m, n, ψ ) converges to 0 for large ψ , e λ U ( m, n, ψ ) still diverges when λ ↑ + ∞ . Onthe other hand, referring to [28] again we have M ( m, n, ψ ) ∼ e ψ ψ m − n Γ( m ) , ψ ↑ + ∞ . Therefore the limit value at x = −∞ does not exist for either e λ M ( m, n, ψ ) or e λ U ( m, n, ψ ). The uniquesolution for the hitting from below case then becomes f ( x ) ≡ . This indicates the LT for the upward FPT does not really exist.On the other hand, consider the solution to FPT from above, i.e. a downward hitting time that theboundary is taken on ∂ D u . By substituting x = + ∞ we get ψ = 0 and λ = −∞ . Refer to [28, Section 13.2 (iii)], depending on the choices of β the limit of U ( m, n, ψ ) at ψ = 0 + hasvarious versions. In order to guarantee the uniqueness of solution we set C = 0. For the limit of M ( m, n, ψ ) we have M ( m, n, ψ ) = 1 + O ( ψ ) , ψ ↓ . The parameter β is involved in the LT. We want a function of solution f ( x ) that is unique in functional forms to all β ≥ ∞ boundary gives the solution f ( x ) = C e λ M ( m, n, ψ ) . Consider the boundary condition on x = a . We write ˆ ψ := (cid:15)e − αa α , ˆ λ := a ( (cid:15)c − (cid:112) c (cid:15) + 2 β ) . Then f ( a ) = 1 gives f ( x ) = e λ M ( m, n, ψ ) e ˆ λ M (cid:16) m, n, ˆ ψ (cid:17) . (27) Remark 4.3
As indicated by our analysis, only a downward LT for the present problem exists. In practicewe are more interested in the burst time of an economic bubble rather than predicting how record-highwould the bubble visit. Therefore the missing solution in upward LT would be a minor issue.
Equation (27) shows the LT for the downward first hitting time. Due to the special function it isdifficult to find the explicit inverse transform. For numerical inversion schemes we refer to [1], wherethree efficient algorithms are provided. However, considering the complicated functional form, it can beimagined that the speed and accuracy in the numerical inverse may not be desired.
The earliest and most successful application of the perturbation technique could be traced back to infinding the solutions of the Schrodinger equation for Hamiltonians of even moderate complexity [37, 42].In mathematical finance, perturbation theory has been studied extensively; see [17, 14, 15]. Inspired byA. Dassios and S. Wu [14], J. Fouque et al. [17], we apply perturbations on the mean-reversion parameter (cid:15) and find the closed-form density for the downward FPT of { X t } t ≥ . W.l.o.g. we let a = 0 and the FPT problem is defined on D u . Consider the function f ∈ C . Assumethere exists a sequence of C functions { f i } i ≥ , such that f = ∞ (cid:88) i =0 (cid:15) i f i . (28)Substitute (28) into (23): ∞ (cid:88) i =0 (cid:15) i A f i = ∞ (cid:88) i =0 (cid:15) i βf i . (29)For any f ∈ C , introduce G f ( x ) := 12 f (cid:48)(cid:48) ( x ) . Rearranging the terms in (29) we have G f − βf + ∞ (cid:88) i =1 (cid:15) i (cid:16) G f i − βf i + ( e − αx − c ) f (cid:48) i − (cid:17) = 0 .
14y assigning proper boundary conditions we split the original Dirichlet problem (23) and (24) intorecursive representations: o (1) : G f − βf = 0 , f ( ∂D ) = (1 , T , (30)and for i ≥ o ( (cid:15) i ) : G f i − βf i + ( e − αx − c ) f (cid:48) i − = 0 , f i ( ∂D ) = (0 , T . (31) Remark 4.4
The o (1) problem in fact is the corresponding boundary value problem for the downwardFPT of Brownian motion. Introduce τ W := inf { t ≥ W t = 0 | W t = x > } . Then f ( x ) = E x (cid:2) e − βτ W (cid:3) . In addition, for i ≥ the function f i embeds the following representation [34] f i ( x ) = E x (cid:20)(cid:90) τ W e − βs (cid:0) e − αW s − c (cid:1) f (cid:48) i − ( W s ) ds (cid:21) . According to [34] the solutions for o (1) and o ( (cid:15) i ) , i ≥ , are unique. Therefore the existence of { f i } i ≥ isguaranteed and the perturbation representation (28) is valid. In the present paper we solve the recursive system up to i = 1 and provide the o ( (cid:15) )-accurate FPTDestimation. Referring to [5], for o (1) we have f ( x ) = e − γx , (32)where γ := √ β . Further let f = f g . Then solving o ( (cid:15) ) we get g ( x ) = γ (cid:0) e − αx − (cid:1) α ( γ + α ) + cx. (33) Proposition 4.5
Let τ ∗ to be the first order approximation of τ . Then the FPTD of τ ∗ is given by P x ( τ ∗ ∈ dt ) = (cid:32) (cid:15) (cid:32) cx + (cid:0) − e − αx (cid:1) ( αt − x )2 αx (cid:33)(cid:33) p ( t ) − (cid:15) α (cid:0) − e − αx (cid:1) e αx ( αt x +1 ) Erfc (cid:32) x √ t + α (cid:114) t (cid:33) , where P x ( · ) = P ( ·| X = x ) = P ( ·|F ) and p ( t ) is the downward FPTD for Brownian motion p ( t ) = x √ π t − e − x t . Erfc ( · ) is the complementary error function given byErfc ( z ) = 2 √ π (cid:90) ∞ z e − y dy. Proof:
Refer to (28), (32) and (33). The first order perturbed LT of the downward FPT is given by f ( x ) = f ( x ) (1 + (cid:15)cx ) + (cid:15)f ( x ) γ (cid:0) e − αx − (cid:1) α ( α + γ ) . (34)15ote that γ = √ β . Therefore f ( x ) is a function of β as well. To emphasize the transform parameterwe denote by f ( β ) and f ( β ) respectively. As mentioned by Remark 4.4, f ( β ) is nothing but the LTfor the downward FPT of Brownian motion. According to [4] this gives p ( t ) := L − { f ( β ) } ( t ) = x √ π t − e − x t . Now consider the inverse transform for the second term in (34). Define˜ l ( β ) := √ βe −√ β αx + √ β . (35)Then the second term can be re-written as f ( β ) γ (cid:0) e − αx − (cid:1) α ( α + γ ) = e − αx − α · ˜ l (2 x β ) . (36)Refer to [4]. The inverse of (35) is given by L − (cid:110) ˜ l ( β ) (cid:111) ( t ) = α x e αx ( αxt +1) Erfc (cid:18) √ t + αx √ t (cid:19) − αxt − √ π t − e − t . (37)According to the property of inverse Laplace transform [31], for constant c L − (cid:26) c ˜ l (cid:18) βc (cid:19)(cid:27) ( t ) = L − (cid:110) ˜ l ( β ) (cid:111) ( ct ) . So let c = x we have L − (cid:110) ˜ l (2 x β ) (cid:111) ( t ) = 12 x L − (cid:110) ˜ l ( β ) (cid:111) (cid:18) t x (cid:19) . (38)Summarizing (36), (37), (38) gives the inverse transform for the second term in (34). This concludes theproof. (cid:3) Remark 4.6
As an approximation, the first order perturbation provides a continuous function but notnecessarily a valid probability density function. In fact, E x [ τ ∗ ] = lim β ↓ f ( β ) = 1 + (cid:15)cx. In the case c > the first order perturbation would provide an extra tiny probability by (cid:15)cx . We willdiscuss the accuracy issue in the later proposition. Follow directly with Proposition 4.5. The tail asymptotics and probability distribution of runningminimum are given explicitly.
Corollary 4.7
The tail asymptotics for P x ( τ ∗ ∈ dt ) are given by P x ( τ ∗ ∈ dt ) ∼ (cid:18) (cid:15) (cid:18) cx − − e − αx α (cid:19)(cid:19) p ( t ) ∼ p ( t ) , t ↓ + , (Left Tail Asymptotics) and P x ( τ ∗ ∈ dt ) ∼ (cid:32) (cid:15) (cid:32) cx + (1 − αx ) (cid:0) − e − αx (cid:1) α x (cid:33)(cid:33) p ( t ) ∼ p ( t ) , t ↑ + ∞ . (Right Tail Asymptotics)16 roof: The left tail asymptotic is given by calculations referring to [9]Erfc (cid:32) x √ t + α (cid:114) t (cid:33) ∼ exp − (cid:32) x √ t + α (cid:114) t (cid:33) , t ↓ + . (39)Consider the right tail. Note that if we repeat using (39), the second term of P x ( τ ∗ ∈ dt ) will remain asa constant while the first term vanishes. This leads to a constant tail asymptotic for t ↑ + ∞ , however, P x ( τ ∗ ∈ dt ) ↓ y ) ∼ e − y y √ π (cid:18) − y (cid:19) , y ↑ + ∞ . (40)Also note the first term of P x ( τ ∗ ∈ dt ) can be re-expressed as (cid:32) (cid:15) (cid:32) cx + (cid:0) − e − αx (cid:1) ( αt − x )2 αx (cid:33)(cid:33) p ( t ) = (cid:32) (cid:15) (cid:32) cx − (cid:0) − e − αx (cid:1) α (cid:33)(cid:33) p ( t ) + (cid:15) (cid:0) − e − αx (cid:1) t x p ( t ) . (41)Then substituting (40) and (41) into P x ( τ ∈ dt ), we find as t ↑ + ∞ , P x ( τ ∗ ∈ dt ) ∼ (cid:32) (cid:15) (cid:32) cx − (cid:0) − e − αx (cid:1) α (cid:33)(cid:33) p ( t ) + (cid:15) (cid:0) − e − αx (cid:1) e − x t √ πt − (cid:15) α (cid:0) − e − αx (cid:1) · √ e − x t √ πtα + (cid:15) α (cid:0) − e − αx (cid:1) · √ e − x t √ πα t = (cid:32) (cid:15) (cid:32) cx − (cid:0) − e − αx (cid:1) α (cid:33)(cid:33) p ( t ) + 0 + (cid:15) − e − αx xα p ( t ) . This completes the proof. (cid:3)
Remark 4.8
Corollary 4.7 indicates the FPTD of τ ∗ has the same tails as the FPTD of Brownianmotion. We know Brownian motion is a null-recurrent Markov process. Therefore we may infer E x [ τ ∗ ] =+ ∞ . Indeed, according to the first moment rule and by (34) we can check E x [ τ ∗ ] = − ∂f ( β ) ∂β (cid:12)(cid:12)(cid:12)(cid:12) β =0 = + ∞ . Corollary 4.9
For fixed t ≥ and a < x , denote the running minimum by X ∗ t := min ≤ u ≤ t { X u } . Also write P x ( τ ∗ ∈ dt ) with parameters (cid:15), x, α, c as p ( t | (cid:15), x, α, c ) . Then the first order perturbed distribution of X ∗ t is given by P x ( X ∗ t ≤ a ) = (cid:90) t p (cid:0) u | (cid:15)e − αa , x − a, α, ce αa (cid:1) du. roof: Introduce { Y t } t ≥ such that Y t = X t − a, ∀ t ≥ . The SDE of Y t is given by dY t = e − αa (cid:15) (cid:0) e − αY t − ce αa (cid:1) dt + dW t , Y = y = x − a. (42)Similarly define Y ∗ t := min ≤ u ≤ t { Y u } . Then P x ( X ∗ t ≤ a ) = P y ( Y ∗ t ≤ . (43)On the other hand, let τ to be the FPT of Y t from y to 0. Note the fact that P y ( Y ∗ t ≤
0) = P y ( τ ≤ t ) . (44)So substituting the parameters in (42) into Proposition 4.5, and considering the equivalence between (44)and (43), we prove the result. (cid:3) Now we consider the accuracy of the perturbation estimation. Denote the actual FPTD of τ by P x ( τ ∈ dt ) . Let the absolute error to be denoted by q τ ( t ) := P x ( τ ∈ dt ) − P x ( τ ∗ ∈ dt ) . (45)Then we show q τ ( t ) is o ( (cid:15) )-accurate. Proposition 4.10
For any t > , there exists a constant M > such that | q τ ( t ) | ≤ M (cid:15) . Moreover, the probabilistic representation of q τ ( t ) is given by q τ ( t ) = (cid:15) E x (cid:20)(cid:90) t ∧ τ (cid:0) e − αX u − c (cid:1) η ( t − u, X u ) du (cid:21) , where η ( t, x ) = − α cosh( αx )2 M ( t, x ) + 1 − e − αx √ πα M ( t, x ) + e − αx √ π M ( t, x ) + c (cid:18) − x t (cid:19) p ( t ) with M ( t, x ) = Erfc (cid:16) x √ t + α (cid:113) t (cid:17) e α t M ( t, x ) = e − x t (cid:2) α t − ( αx + 1) t + x (cid:3) t − M ( t, x ) = e − x t ( αt − x ) t − . roof: To emphasize the dual effects of f as a function of x, β , denote by f ( β, x ) := f ( x ) = f ( β ).Let η ( t, x ) := L − (cid:8) ∂∂x f ( β, x ) (cid:9) ( t ) to be the inverse transform of f (cid:48) ( x ). Consider using same tricks inthe proof of Proposition 4.5. After standard calculations we show η ( t, x ) embeds the explicit form as inabove.Now we prove the probabilistic representation and uniform boundness. Let h ( x ) := ( e − αx − c ). Bythe solution of η ( t, x ) we know lim t ↑ + ∞ | h ( x ) η ( t, x ) | = 0 , ∀ x ∈ D u . On the other hand, according to Corollary 3.11, P x ( τ < + ∞ ) = 1. Thereforelim t ↑ + ∞ (cid:90) t ∧ τ | h ( X u ) η ( t − u, X u ) | du = (cid:90) τ lim t ↑ + ∞ | h ( X u ) η ( t − u, X u ) | du = 0 . Since (cid:82) t ∧ τ | h ( X u ) η ( t − u, X u ) | du is continuous on t andlim t ↓ (cid:90) t ∧ τ | h ( X u ) η ( t − u, X u ) | du = 0 , so there exists M > (cid:90) t ∧ τ | h ( X u ) η ( t − u, X u ) | du ≤ M, ∀ t ≥ , and { X u } ≤ u ≤ t ∈ D u . (46)Let ˜ q τ ( t ) defined by ˜ q τ ( t ) = (cid:15) E x (cid:20)(cid:90) t ∧ τ h ( X u ) η ( t − u, X u ) du (cid:21) . By (46) we immediately have | ˜ q τ ( t ) | ≤ M (cid:15) . (47)For β ∈ C and Real( β ) ≥
0, consider the Laplace transform of ˜ q τ ( t ) L { ˜ q τ ( t ) } ( β ) = (cid:15) (cid:90) ∞ e − βt E x (cid:20)(cid:90) t ∧ τ h ( X u ) η ( t − u, X u ) du (cid:21) dt. Based on (46) and the dominated convergence theorem, we change the order of integral and expectation
L { ˜ q τ ( t ) } ( β ) = (cid:15) E x (cid:20)(cid:90) ∞ (cid:90) t ∧ τ e − βt h ( X u ) η ( t − u, X u ) dudt (cid:21) . (48)In addition, (46) also gives the Fubini’s theorem so (cid:90) ∞ (cid:90) t ∧ τ e − βt h ( X u ) η ( t − u, X u ) dudt = (cid:90) ∞ (cid:90) τ { u ≤ t } e − βt h ( X u ) η ( t − u, X u ) dudt = (cid:90) τ (cid:90) ∞ { u ≤ t } e − βt h ( X u ) η ( t − u, X u ) dtdu = (cid:90) τ L (cid:8) { u ≤ t } η ( t − u, X u ) (cid:9) ( β ) h ( X u ) du. Note { u ≤ t } is the indicator function, which can also be written as the Heaviside step function H u ( t ).Consider the fact [31] that L { H u ( t ) η ( t − u, x ) } ( β ) = e − βu L { η ( t, x ) } ( β ) , L { η ( t, x ) } ( β ) = L (cid:26) L − (cid:26) ∂∂x f ( β, x ) (cid:27) ( t ) (cid:27) ( β ) = ∂∂x f ( β, x ) . Therefore (48) can be re-expressed as
L { ˜ q τ ( t ) } ( β ) = (cid:15) E x (cid:20)(cid:90) τ e − βu h ( X u ) ∂∂x f ( β, X u ) du (cid:21) . (49)In the next step we show L { ˜ q τ ( t ) } ( β ) indeed is the LT for the error function q τ ( t ). Then theuniqueness of inverse LT concludes our proof. To see this, let Q ( β, x ) := f ( β, x ) − f ( β, x ) , (50)where follow similar convention f ( β, x ) := f ( x ) and f ( x ) is as introduced in (34). Note that f ( β, x ) = L { P x { τ ∈ dt }} ( β ) and f ( β, x ) = L (cid:8) P x { τ ∗ ∈ dt } (cid:9) ( β ). By the linearity of LT we therefore have Q ( β, x ) = L { q τ ( t ) } ( β ) . (51)As f, f ∈ C , so is Q ( β, x ) ∈ C . Apply the infinitesimal generator A on Q ( β, x ). Then by (23), (30)and (31) with i = 1, after standard calculations we get A Q − βQ = − (cid:15) hf (cid:48) , x ∈ D u . (52)Note (52) is an equation about x and f (cid:48) is a short for ∂∂x f ( β, x ). Since f and f share the same boundaryconditions, so the boundary condition of ODE (52) is given by Q ( ∂ D u ) = (0 , T . (53)According to [34], the boundary value problem (52) and (53) has the following unique solution Q ( β, x ) = (cid:15) E x (cid:20)(cid:90) τ e − βu h ( X u ) ∂∂x f ( β, X u ) du (cid:21) . The uniqueness in ODE solution and the uniqueness in inverse LT indicates q τ ( t ) = ˜ q τ ( t ) . (cid:3) Remark 4.11
Proposition 4.10 shows the error bound is uniformly valid on t ≥ . When (cid:15) ↓ + , theerror converges to . This is true as when X t → W t pathwisely, referring to Proposition 4.5 we have P x ( τ ∗ ∈ dt ) → p ( t ) , ∀ t ≥ . Remark 4.12
On the other hand, the conclusion in Proposition 4.10 does not restrict applying perturba-tion for (cid:15) > . In fact, an exact error function is given and we can estimate the error level via simulation.Even in the case that (cid:15) > , there are possibilities that (cid:12)(cid:12)(cid:12)(cid:12) E x (cid:20)(cid:90) t ∧ τ (cid:0) e − αX u − c (cid:1) η ( t − u, X u ) du (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) << (cid:15) , t ∈ (0 , + ∞ ) . Model Implementation
For practical purpose it is more interesting to take the volatility into account. We extend SDE (1) byadding a constant volatility σ > dX t = (cid:15) ( e − αX t − c ) dt + σdW t , X = x ∈ R . (54)Introduce the scaled version of { X t } t ≥ and define (cid:110) ˜ X t (cid:111) t ≥ as˜ X t := X t σ . Then by setting ˜ (cid:15) := (cid:15)σ , ˜ α := ασ and ˜ x = xσ , we see (cid:110) ˜ X (cid:111) t ≥ indeed is the diffusion process described bySDE (1): d ˜ X t = ˜ (cid:15) (cid:16) e − α ˜ X t − c (cid:17) dt + dW t , ˜ X = ˜ x ∈ R . In addition, for a ≤ x , by letting ˜ a = aσ the result of running minimum in Corollary 4.9 can be extendedaccordingly. In this section we provide a calibration scheme for the extended SDE (54). Denote the observations ofasset prices { P t } t =0 , ,...,N by P t = P e ˆ X t , t = 0 , , ..., N. (55)Then (cid:110) ˆ X t (cid:111) t =0 ,...,N represents the normalised log-price with ˆ X = 0. Let { ˆ r t } t =1 ,...,N to be the log-returnof { P t } t =0 , ,...,N . By definition we haveˆ r t = ˆ X t − ˆ X t − , t = 1 , ..., N. (56)Consider the calibration based on { ˆ r t } t =1 ,...,N . Mathematically, there are 4 parameters to be decided.Therefore at least 4 different statistical quantities should be provided. A natural candidate is the first fourmoments of { ˆ r t } t =1 ,...,N . However, on the one hand, as we discussed in Section 2, the bubble dynamicin different regimes could have totally different statistical behaviours. So global moments on the wholetime-series may not be representative. On the other hand, from Proposition 3.4, { X t } t ≥ has a verycomplicated probability density. Following the proposition, we cannot easily get the explicit expressioneven for the first moment. Instead of using traditional moments calibration, we provide an alternativescheme with the piecewise time-series under different bubble regimes.Recall those three regimes of { X t } t ≥ in a bubble cycle, according to which we make the followingassumptions: • Regime I), displacement. During this period we assume X t ≈
0. SDE (54) then can be simplifiedas dX t ≈ (cid:15) (1 − c ) dt + σdW t . (57)21 Regime II), boom. In this stage the dynamic follows SDE (54) but will visit the equilibrium level.Denote the level by X R , we have e − αX R = c. (58) • Regime III) euphoria (& profit taking). Within these two steps X t hits the record-high level.Assume e − αX t ≈ dX t ≈ − c(cid:15)dt + σdW t . (59)Besides, we further assume each regime could be recognised from the data. Let { , , ..., t } , { t , ..., t } , { t , ..., t } to be the time periods for regimes I, II and III. Then denote the piecewise time-series in eachregime by ˆ X I := (cid:110) ˆ X t (cid:111) t =0 ,...,t , ˆ X II := (cid:110) ˆ X t (cid:111) t = t ,...,t , ˆ X III := (cid:110) ˆ X t (cid:111) t = t ,...,t . The corresponding time-series for log-returns are given byˆ r I := { ˆ r t } t =1 ,...,t , ˆ r II := { ˆ r t } t = t +1 ,...,t , ˆ r III := { ˆ r t } t = t +1 ,...,t . Also assume that the equilibrium level is observable and denote the observation byˆ X R . We now consider parameter estimates. Start with ˆ (cid:15) and ˆ c . The general idea is to take the expectedlog-returns from regimes I and III in to account. Let¯ r I := Mean (cid:0) ˆ r I (cid:1) and ¯ r III := Mean (cid:0) ˆ r III (cid:1) to be the annualized sample averages of returns. By matching the sample means with theoretical expec-tations from dX t in (57) and (59), we have the following equations ˆ (cid:15) (1 − ˆ c ) = ¯ r I − ˆ (cid:15) ˆ c = ¯ r III . Solving the equations we get ˆ (cid:15) = ¯ r I − ¯ r III ˆ c = − ¯ r III ¯ r I − ¯ r III . (60) Remark 5.1
Note that according to our assumptions, regime I should provide positive trend ( ¯ r I ≥ )while regime III generates negative moves ( ¯ r III ≤ ). Therefore ˆ (cid:15) and ˆ c are guaranteed to be positive.Moreover, since ≤ − ¯ r III ≤ ¯ r I − ¯ r III , so ≤ ˆ c ≤ . Remark 5.2
In order to have a more effective calibration, in ¯ r I and ¯ r III estimations we can (*) takethe average of only positive returns in regime I and only the negative returns in regime III. In addition,we are more interested in the longer term trend rather than the daily trend. So (**) using monthly rollingreturns would help in enhancing the estimation stability. We add (*) and (**) as special data cleaningtreatments in our algorithm. σ . By observing (57) and (59) we see the volatilities in ˆ r I and ˆ r III are provided by theBrownian motion part only. Let ˆ r I & III := ˆ r I ∪ ˆ r III . Then we can compute ˆ σ byˆ σ = StdDev (cid:0) ˆ r I & III (cid:1) . As an alternative plan, notice that usually regime III has more volatile time-series. Therefore in order tocapture a more significant volatility we choose to use ˆ r III only:ˆ σ = StdDev (cid:0) ˆ r III (cid:1) . (61)Given ˆ c , the last parameter ˆ α is easy to compute. Based on (58), we immediately haveˆ α = − ln (ˆ c )2 ˆ X R . (62)We summarise the calibration algorithm in Algorithm 1. Algorithm 1 { X t } t ≥ Parameter Calibration1. Determine the time ranges for regimes I-III, and correspondingly calculate the floored log-priceˆ X I , ˆ X II , ˆ X III by (55). Identify the equilibrium level ˆ X R .2. Calculate the monthly rolling log-returns of ˆ r Im and ˆ r IIIm from ˆ X I and ˆ X III respectively. Calculate¯ r I and ¯ r III via ¯ r I = M ean (cid:18) ˆ r Im (cid:12)(cid:12)(cid:12)(cid:12) ˆ r Im ≥ (cid:19) × r III = M ean (cid:18) ˆ r IIIm (cid:12)(cid:12)(cid:12)(cid:12) ˆ r IIIm ≤ (cid:19) × , and use Equation (60) to calibrate ˆ (cid:15), ˆ c .3. Calculate the daily log-return time-series ˆ r IIId from ˆ X III . Compute annualised return ˆ r III viaˆ r III = ˆ r IIId × √ , and calibrate ˆ σ using Equation (61).4. Substitute ˆ X R from step 1 and ˆ c from step 2 into equation (62) to calibrate ˆ α . Remark 5.3
We need to highlight that Algorithm 1 relies on two judgmental decisions, i.e. 1) thetime range for different regimes and 2) the equilibrium level. During an asset price increasing period(displacement, boom, euphoria), from the economical point of view it is not difficult to differentiate thosethree regimes. Even when there is no clear economical signal, we can still split the time-series equallyinto three pieces. However, for ˆ X R , without a significant price drop, mathematically it is very challengingto decide the equilibrium level. Therefore fundamental analysis from economics may be required. Theenhancement of Algorithm 1 will be remained in the future work. .3 Numerical Examples We provide three numerical examples. The first two exercises are similar in nature, where based onhistorical data we verified the effectiveness of { X t } t ≥ in capturing bubble dynamics. In the third exercisewe predicted drop-down probabilities for BitCoin. The US dot-com bubble [8] could be observed from the technology-dominated NASDAQ Composite Index(US ticker symbol ˆIXIC). From mid 90’s, ˆIXIC grew exponentially from below 1,000 USD to about5,000 USD. The index hit its historical maximum in 2000-03-10, and at which date the total tradingamount exceeded 10 Trillion USD (according to Yahoo Finance). After then the market collapsed rapidlyand dropped back to about 1,000 USD in 2002.In this exercise we used the adjusted daily close price of ˆIXIC from 1997-01-02 to 2003-12-30. Thedata was downloaded from Yahoo Finance. Note that, for the purpose of burst time prediction, there isno sense to calibrate the model using the full-cycle data. Therefore only a truncated time-series was usedin model calibration. To be more specific, we chose the calibration regimes as follows ˆ X I : 1997-01-02 ( P = 1 , P t = 1 , X II : 1997-06-26 ( P t = 1 , P t = 2 , X III : 1999-02-10 ( P t = 2 , P t = 3 , (cid:110) ˆ X t (cid:111) t =0 , ,...,N . By observation we set the equilibriumlevel to be ˆ X R = 0 .
67 ( P R = 2 , { X t } t ≥ calibration, i.e. from 1997-01-02 to 2000-10-18. For the MLE OU calibration algorithm, cf.[41]. We estimated the mean and volatility directly in the DBM. 1,000 paths between 1997-01-02 and2003-12-30 were simulated by three different models. In Figure 3, apart from the historical price of ˆIXIC,we demonstrate the best path among 1,000 simulations for each model. It is clear by the graph that ournew model provides better fit than existing models. To quantitatively see the closeness of different pathsto the historical dynamic, the correlations for each model were calculated: { X t } t ≥ : 91 . , OU : 81 . , DBM : 72 . . As expected, { X t } t ≥ provided the highest correlation while DBM was the worst among three models.To further explain our algorithm, we plot calibration regimes in Figure 4. We also show 10,000 onwardsimulation paths for { X t } t ≥ with X = ˆ X t . From the figure we see the historical prices are fully coveredby the simulation paths. This indicates our model is effective.24igure 3: Model calibration comparisons for NASDAQ index (US ticker symbol ˆIXIC) from 1997-01-02to 2003-12-30. Red curve: historical adjusted log-price; blue curve: the best of 1,000 simulations from { X t } t ≥ ; orange curve: the best of 1,000 simulations from OU process; green curve: the best of 1,000simulations from DBM. Calibration parameters, { X t } t ≥ : (ˆ (cid:15), ˆ α, ˆ σ, ˆ c ) = (0 . , . , . , . OU :(ˆ κ, ˆ µ, ˆ σ ) = (0 . , . , . BM : (ˆ µ, ˆ σ ) = (0 . , . X = ˆ X t .Green zone indicates regime I, the displacement stage; yellow zone indicates regime II, the boom stage;red zone indicates regime III, the euphoria & profit taking stages. Blue curve shows the historical dataused for calibration. Red curve, covered by shadowed region, shows the historical data after t . Theshadowed region plots 10,000 simulation paths. 25 .3.2 2006-01-04 to 2008-12-31 Shanghai Stock Exchange Composite Index We use a second example to confirm our observations from Section 5.3.1. The 2007 Chinese stock marketcrash [24] just happened before the 2008 global financial crisis. Starting in early 2006, the ShanghaiStock Exchange Composite Index (US ticker symbol SSEC, China ticker symbol 000001.SS) increasedfrom about 1,000 CNY to 6,092 CNY in mid-October, 2007. And within one year’s time, from October2007 to October 2008, the price dropped below 1,800 CNY. Similar to the pattern in ˆIXIC, the historicallog-price of SSEC dropped rapidly after the sudden peak, and before which there was a sharp increase.The exercise settings were the same as in Section 5.3.1. We only mention the regime settings andmake comments where are necessary. ˆ X I : 2006-01-04 ( P = 1 , P t = 1 , X II : 2006-03-06 ( P t = 1 , P t = 4 , X R = 1 .
23 ( P R = 4 , X III : 2007-05-30 ( P t = 4 , P t = 3 , { X t } t ≥ : 96 . , OU : 88 . , DBM : 84 . . Similar plot for the algorithm illustration and 10,000 simulation paths is given in Figure 6. Through thisexercise we further confirm that our new model is a good candidate for describing economic bubbles.26igure 5: Model calibration comparisons for Shanghai Stock Exchange Composite index (US ticker symbolSSEC, China ticker symbol 000001.SS) from 2006-01-04 to 2008-12-31. Red curve: historical adjusted log-price; blue curve: the best of 1,000 simulations from { X t } t ≥ ; orange curve: the best of 1,000 simulationsfrom OU process; green curve: the best of 1,000 simulations from DBM. Calibration parameters, { X t } t ≥ :(ˆ (cid:15), ˆ α, ˆ σ, ˆ c ) = (0 . , . , . , . OU : (ˆ κ, ˆ µ, ˆ σ ) = (3 . , . , . BM : (ˆ µ, ˆ σ ) = (0 . , . X =ˆ X t . Green zone indicates regime I, the displacement stage; yellow zone indicates regime II, the boomstage; red zone indicates regime III, the euphoria & profit taking stages. Blue curve shows the historicaldata used for calibration. Red curve, covered by shadowed region, shows the following historical dataafter t . The shadowed region plots 10,000 simulation paths.27 .3.3 BitCoin Downward Probability Estimation ˆ X I : 2016-01-01 ( P = 433) to 2016-05-30 ( P t = 528);ˆ X II : 2016-05-30 ( P t = 528) to 2017-08-13 ( P t = 4 , X R = 2 .
30 ( P R = 4 , X III : 2017-08-13 ( P t = 4 , P t = 14 , (cid:15) = 0 .
51; ˆ α = 0 .
08; ˆ σ = 0 .
91; ˆ c = 0 . . (66)The prediction was made on 2017-12-10 with the price at P t = 14 , . We considered 0% to 60%drops from P t . To convert the drop percentages from { P t } t =0 , ,...,N to the log-price space, { X t } t =0 , ,...,N , Note that, the data in our record does not correspond to the close price of 2017-12-10. In fact, the data was downloadedwhen the market was still under trading.
28e calculated the hitting level a via (55). Referring to Section 5.1, we transferred parameters in (66),together with a , from the extended SDE (54) to the parameters in the standard SDE (1). By Corollary4.9, in the end we were able to have the probability distribution for the minimum price within onemonth time. On the other hand, an error evaluation on the perturbed FPTD should be given. Refer toProposition 4.10. The relative error is given by e ( t ) := (cid:12)(cid:12)(cid:12)(cid:12) P ( τ ∗ ∈ dt ) P ( τ ∗ ∈ dt ) + q τ ( t ) (cid:12)(cid:12)(cid:12)(cid:12) . Using the probabilistic representation, we estimated q τ ( t ) via 10,000 paths simulation. It should benoticed that, the relative error generally is high at tails as the actual density converges to 0. Therefore itis not necessary to compute the relative error at each point. In fact, we are more concerned that whetherthe peak of the distribution would be changed by perturbations. So only relative errors on the densitypeak were computed. Table 1 summarises the results.Percentage of Drop Price P l (USD) Probability P ( P ∗ t ≤ P l ) Peak Relative Error0% 14,371.62 100.00% 0.00%5% 13,653.05 84.85% 4.97%10% 12,934.47 69.38% 1.48%15% 12,215.89 54.25% 0.90%20% 11,497.30 40.19% 0.04%25% 10,778.72 27.88% 0.59%30% 10,060.14 17.87% 0.86%35% 9,341.56 10.38% 0.88%40% 8,622.98 5.35% 0.68%45% 7,904.40 2.37% 1.69%50% 7,185.81 0.86% 1.45%55% 6,467.23 0.25% 2.40%60% 5,748.65 0.05% 1.78%Table 1: BitCoin downward price prediction between 2017-12-10 and 2018-01-12. Columns 1-4 correspondto the percentages of price drop, dropped price P l , probability of the lowest price P ∗ t below P l , and therelative error in density peaks.First by checking the last column (relative errors), we see in general the perturbation model is accurate.The largest error was in the 5% drop. In this case the hitting level is very close to the initial price P t .As a result, the density curve will shrink to the y-axis. Therefore a larger error is expected. Analogously,large errors might also exist in the case that hitting levels are far to the initial price. By ruling out theextreme drops, in the range of 10% to 50%, we see the estimation errors remained below 2%.We now consider the possibility of market collapse. Referring to the scenarios in ˆIXIC and 000001.SS,we found their largest drops in a month were about 30%, and which happened in the spring of 2000 andautumn of 2008, respectively. Then check the probability of 30% drop for BitCoin. From Table 1 we onlysee about 17.87%. In fact, even for a 20% drop, the probability was about 40.19%. This means there29as more than half chance that the price would remain above 11,497.30. Therefore we concluded thatthe market was unlikely to collapse in the next month.We collected the one month data from 2017-12-10 to 2018-01-12 and plot the time-series in Figure 8.From the graph we see the lowest close price was 12,531.52 on 2017-12-30. This verified our conclusionthat the market would not collapse immediately. On the other hand, compare the probabilities in Table1 with the thresholds in Figure 8. The price on 2017-12-30 broke the 10% drop threshold, where theprobability given by our prediction was 69.38%. This further confirmed that the model is effective.Figure 8: BitCoin close price between 2017-12-10 and 2018-01-12. The probabilities of different thresholdsare reported in Table 1. In this paper we find a new diffusion process which can be used in modelling economic bubbles. Thesimple form of the model enables us to deduce its downward FPTD explicitly. Therefore the paperprovides a useful tool in estimating the burst time of an economic bubble. Numerical examples in Section5 consistently confirm that the model and its prediction are effective. Results in Section 3 show theprocess has desirable properties which potentially can be employed in the future option pricing work.The perturbation technique, as introduced in Section 4, can be extended in finding explicit FPTDs ofother diffusion processes. In another working paper of us, the corresponding closed-form FPTDs havebeen found for the OU and Bessel processes. One remaining issue is the exact simulation for the process.This requires further understandings to the θ ( r, s ) function and we leave it for the future work.30 eferences [1] Joseph Abate and Ward Whitt. A unified framework for numerically inverting Laplace transforms. INFORMS Journal on Computing , 18(4):408–421, 2006.[2] Milton Abramowitz and Irene A Stegun.
Handbook of mathematical functions: with formulas, graphs,and mathematical tables , volume 55. Courier Corporation, 1964.[3] Pauline Barrieu, A Rouault, and M Yor. A study of the HartmanWatson distribution motivatedby numerical problems related to the pricing of Asian options.
Journal of Applied Probability ,41(4):1049–1058, 2004.[4] Harry Bateman.
Tables of integral transforms [volumes I & II] , volume 1. McGraw-Hill BookCompany, 1954.[5] Andrei N Borodin and Paavo Salminen.
Handbook of Brownian motion-facts and formulae .Birkh¨auser, 2012.[6] Chris Brooks and Apostolos Katsaris. A three-regime model of speculative behaviour: modelling theevolution of the S&P 500 composite index.
The Economic Journal , 115(505):767–797, 2005.[7] Chris Brooks and Apostolos Katsaris. Trading rules from forecasting the collapse of speculativebubbles for the S&P 500 composite index.
The Journal of Business , 78(5):2003–2036, 2005.[8] John Cassidy. dotcom: how America lost its mind and money in the internet era, 2003.[9] Marco Chiani, Davide Dardari, and Marvin K Simon. New exponential bounds and approxima-tions for the computation of error probability in fading channels.
IEEE Transactions on WirelessCommunications , 2(4):840–845, 2003.[10] Alexander MG Cox and David G Hobson. Local martingales, bubbles and option prices.
Financeand Stochastics , 9(4):477–492, 2005.[11] AB Olde Daalhuis. Confluent hypergeometric functions.
NIST Handbook of Mathematical Functions,FWJ Olver, DW Lozier, RF Boisvert, and CW Clark, eds., Cambridge University, New York , pages321–349, 2010.[12] Mike Dash.
Tulipomania: The story of the world’s most coveted flower and the extraordinary passionsit aroused . Hachette UK, 2011.[13] Angelos Dassios and Jayalaxshmi Nagaradjasarma. The square-root process and Asian options.
Quantitative Finance , 6(4):337–347, 2006.[14] Angelos Dassios and Shanle Wu. Perturbed Brownian motion and its application to Parisian optionpricing.
Finance and Stochastics , 14(3):473–494, 2010.[15] Peter W Duck, Chao Yang, David P Newton, and Martin Widdicks. Singular perturbation techniquesapplied to multiasset option pricing.
Mathematical Finance , 19(3):457–486, 2009.3116] Vladimir Filimonov, Guilherme Demos, and Didier Sornette. Modified profile likelihood inferenceand interval forecast of the burst of financial bubbles.
Quantitative finance , 17(8):1167–1186, 2017.[17] Jean-Pierre Fouque, George Papanicolaou, Ronnie Sircar, and Knut Sølna.
Multiscale stochasticvolatility for equity, interest rate, and credit derivatives . Cambridge University Press, 2011.[18] Peter M Garber.
Famous first bubbles: The fundamentals of early manias . mit Press, 2001.[19] Philip Hartman and Geoffrey S Watson. ” Normal” distribution functions on spheres and the modifiedBessel functions.
The Annals of Probability , pages 593–607, 1974.[20] Steven L Heston, Mark Loewenstein, and Gregory A Willard. Options and bubbles.
The Review ofFinancial Studies , 20(2):359–390, 2006.[21] Hana Hlad´ıkov´a and Jarmila Radov´a. Term structure modelling by using Nelson-Siegel model.
European Financial and Accounting Journal , 7(2):36–55, 2012.[22] Jeff Holt. A summary of the primary causes of the housing bubble and the resulting credit crisis: Anon-technical paper.
The Journal of Business Inquiry , 8(1):120–129, 2009.[23] Kazuyuki Ishiyama. Methods for evaluating density functions of exponential functionals representedas integrals of geometric Brownian motion.
Methodology and Computing in Applied Probability ,7(3):271–283, 2005.[24] Zhi-Qiang Jiang, Wei-Xing Zhou, Didier Sornette, Ryan Woodard, Ken Bastiaensen, and PeterCauwels. Bubble diagnosis and prediction of the 2005–2007 and 2008–2009 Chinese stock marketbubbles.
Journal of economic behavior & organization , 74(3):149–162, 2010.[25] Cassidy John. Dot. con: How america lost its mind and money in the internet era, 2003.[26] Ioannis Karatzas and Steven Shreve.
Brownian motion and stochastic calculus , volume 113. SpringerScience & Business Media, 2012.[27] Alexander Kiselev and Lenya Ryzhik. A simple model for asset price bubble formation and collapse. arXiv preprint arXiv:1009.0299 , 2010.[28] Daniel W Lozier. NIST digital library of mathematical functions.
Annals of Mathematics andArtificial Intelligence , 38(1-3):105–119, 2003.[29] Hiroyuki Matsumoto, Marc Yor, et al. Exponential functionals of Brownian motion, i: Probabilitylaws at fixed time.
Probability surveys , 2:312–347, 2005.[30] Hyman P Minsky and Henry Kaufman.
Stabilizing an unstable economy , volume 1. McGraw-HillNew York, 2008.[31] Fritz Oberhettinger and Larry Badii.
Tables of Laplace transforms . Springer Science & BusinessMedia, 2012. 3232] Frank WJ Olver.
NIST handbook of mathematical functions hardback and CD-ROM . Cambridgeuniversity press, 2010.[33] Goran Peskir. On the fundamental solution of the KolmogorovShiryaev equation. In
From stochasticcalculus to mathematical finance , pages 535–546. Springer, 2006.[34] Goran Peskir and Albert Shiryaev.
Optimal stopping and free-boundary problems . Springer, 2006.[35] Peter CB Phillips, Shuping Shi, and Jun Yu. Testing for multiple bubbles: Limit theory of real-timedetectors.
International Economic Review , 56(4):1079–1134, 2015.[36] Carmen M Reinhart and Kenneth S Rogoff. Is the 2007 US sub-prime financial crisis so different?An international historical comparison.
American Economic Review , 98(2):339–44, 2008.[37] Erwin Schr¨odinger. Quantisierung als eigenwertproblem.
Annalen der physik , 385(13):437–490, 1926.[38] Robert C Shiller. Irrational exuberance.
Philosophy & Public Policy Quarterly , 20(1):18–23, 2000.[39] Albert N Shiryaev. Quickest detection problems in the technical analysis of the financial data. In
Mathematical FinanceBachelier Congress 2000 , pages 487–521. Springer, 2002.[40] AN Shiryaev. The problem of the most rapid detection of a disturbance in a stationary process. In
Soviet Math. Dokl , volume 2, 1961.[41] William Smith. On the simulation and estimation of the mean-reverting Ornstein-Uhlenbeck process.
Commodities Markets and Modelling , 2010.[42] Carlos E Soliverez. General theory of effective hamiltonians.
Physical Review A , 24(1):4, 1981.[43] Denis Talay. Numerical solution of stochastic differential equations. 1994.[44] Christopher Wood.
The bubble economy: the Japanese economic collapse . Sidgwick & Jackson, 1992.[45] Marc Yor. On some exponential functionals of Brownian motion.