Likelihood-free Bayesian inference for alpha-stable models

Lik eliho o d-free Ba y esian inference for α -stable mo dels G. W. Pe ters ∗ † , S.A. Si sson † and Y. F an † June 2, 2 0 18 Abstract α -stable distributions are utilised as mo dels for hea vy-tail ed n oise in man y areas of statistics, finance and signal pr o cessing engineering. Ho we ver, in gen- eral, neither univ ariate n or multiv a r iate α -stable mo dels admit closed f orm densities whic h can b e ev a lu ated p oin twise. This complicates the inferentia l pro cedure. As a resu lt, α -stable mo dels are practically limited to the univ ari- ate setting under the Ba yesian paradigm, and to biv ariate mo dels under the classical framework. In this article w e dev elop a nov el Ba yesia n approac h to mo delling un iv ariate and multiv a r iate α -stable distribu tions based on recen t adv ances in “like liho o d-free” in ference. W e presen t an ev aluati on of the p er- formance of this pro cedu r e in 1, 2 and 3 dimen s ions, and provide an analysis of real daily currency exc hange r ate d ata. The prop osed approac h pro vid es a feasible inferen tial metho d ology at a mod erate computatio n al cost. Keyw ords: α -stable distributions; Appro ximate Ba y esian compu tation; Lik eliho o d- free inference; Sequen tial Mon te Carlo samplers. ∗ Corresp onding Author. Email: Ga rethPeters@unsw.edu.au † School of Mathematics and Statistics, Univ ersity of New South W ales, Sydney 2052, Australia. 1 1 In tro ducti on Mo dels constructed with α -stable distributions p ossess sev eral useful pro p erties, in- cluding infinite v aria nce, sk ewness and heavy tails (Zolotarev 1986; Alder et al. 1998; Samoro dnitsky and T aqqu 1994; Nolan 2007). α -stable distributions pr ovide no gen- eral analytic expressions for the densit y , median, mo de or entrop y , but are uniquely sp ecified b y their c hara cteristic function, whic h has sev eral parameterizations. Con- sidered as generalizations o f the Gaussian distribution, they are defined as the class of lo cation-scale distributions whic h are closed under con v olutions. α - stable distribu- tions ha ve found application in man y areas of statistics, finance and signal pro cessing engineering as mo dels for impulsiv e, hea vy ta iled noise pro cesses (Mandelbrot 1960; F ama 1965; F ama and Roll 196 8; Nikias and Shao 1995; Go dsill 2000; Melc hiori 20 06). The univ ariate α -stable distribution is typ ically sp ecified by four parameters: α ∈ (0 , 2] determining the rate of ta il deca y; β ∈ [ − 1 , 1] determining the degree and sign of asymmetry (sk ewness); γ > 0 the scale (under some parameterizations); and δ ∈ R the lo cation (L evy 1924). The para meter α is termed the c hara cteristic exp onen t, with small and large α implying heav y and ligh t tails resp ectiv ely . G a ussian ( α = 2 , β = 0) and Cauch y ( α = 1 , β = 0) distributions provide the only analytically tractable sub- mem b ers of this family . In general, as α -stable mo dels admit no closed form expres sion for the densit y whic h can b e ev a luated p oin twis e (excepting Gaussian and Cauc h y mem b ers), inference t ypically pro ceeds via the c haracteristic function. This pap er is concerned with constructing b oth univ ariate and m ultiv aria te Bay esian mo dels in whic h the lik eliho o d mo del is from the class of α -stable distributions. This is know n to b e a difficult problem. Existing metho ds for Ba y esian α - stable mo d- els are limited to the univ ariate setting (Buc kle 1995; Go dsill 19 99; Go dsill 2000; Lom bardi 2007; Casarin 2004; Salas-Gonzalez et al. 2006). 2 Inferen tial pro cedures fo r α -stable mo dels ma y b e classified a s auxiliary v ar ia ble metho ds, in ve rsion plus series expansion approaches and densit y estimation metho ds. The auxiliary v ariable Gibbs sampler (Buc kle 1995) increases the dimension of the parameter space fro m 4 ( α, β , γ and δ ) to n + 4, where n is the num b er of observ ations. As strong correlations b etw een para meters and large sample sizes are common in the α -stable setting, this results in a slo wly mixing Mark ov c hain since Gibbs mov es are limited to mo ving parallel to the axes (e.g. Neal 2003). Other Mark ov c hain Mon te Carlo (MCMC) samplers (DuMouc hel 1975; Lom bardi 2007) adopt inv ersion tec hniques for nume r ical in t egr a tion of the c haracteristic function, emplo ying in v erse F ourier tra nsforms com bined with a series expansion (Bergstrom 1953) to accurately estimate distributiona l t a ils. This is p erfo rmed a t eac h iteration of the Mark ov c hain to ev aluate the lik eliho o d, and is accordingly computationally intens ive. In addition the qualit y of the resulting approxim a tion is sensitiv e to the spacing of the fast F ourier transform grid and the p oint at whic h the series expansion b egins (Lom bardi 2007). Univ ariate density estimation metho ds include inte g ral r epresen ta t io ns ( Zolotarev 198 6), parametric mixtures (Nolan 1997; Mc Cullo c h 1998) a nd n umerical estimation through splines and series expansions (Nolan et al. 2001; Nolan 2008). McCullo c h (199 8) a p- pro ximates symmetric stable distributions using a mixture of G aussian and Cauch y densities. Doganoglu and Mittnik (1998) and Mittnik and Rac hev (1991) approxi- mate the α -stable densit y through spline p olynomials, and Kuruoglu et al. (1997) via a mixture of G aussian distributions. P arameter estimation has b een p erfo rmed b y an exp ectation-maximization (EM) a lg orithm (Lom bardi and G o dsill 2 006) and b y metho d of (par t ial) momen ts (Press 1972; W eron 2006). Implemen ted within an MCMC sampler, suc h densit y estimation metho ds would b e highly computat io nal. None of the a b o ve metho ds eas ily generalize to the m ultiv ariate setting. It is curren tly only practical to numeric a lly ev alua t e t wo dimensional α - stable densities via 3 in v ersion of the c haracteristic function. Here the required computation is a f unction of α and t he num b er and spread of masses in the discrete sp ectral represen ta tion (Nolan et al. 2001; Nolan 1997). Bey ond tw o dimensions t his pro cedure b ecomes un tenably slo w with limited a ccuracy . In this article w e dev elop practical Ba ye sian inferen tial metho ds to fit univ a riate and m ultiv ar iate α - stable mo dels. T o the b est of our kno wledge, no pra ctical Bay esian metho ds hav e b een dev elop ed for the m ultiv ariate mo del as the required computa- tional complexit y increases dramatically with mo del dimension. The same is true of classical metho ds b ey ond t wo dimensions. Our approac h is based on recen t dev elop- men ts in “lik eliho o d-free” inference, whic h p ermits appro ximate po sterior sim ulation for Bay esian mo dels without the need to explicitly ev aluate the lik eliho o d. In Section 2 w e briefly in tro duce like liho o d-free inference and the sampling frame- w ork used in this art icle. Section 3 presen ts the Bay esian α - stable mo del, with a par- ticular fo cus on summary statistic sp ecification, a critical comp onent of like liho o d-fr ee inference. W e provide an ev aluation of the p erfo r ma nce of the pro p osed metho dology in Section 4, based on controlled sim ulatio n studies in 1, 2 and 3 dimensions. Fi- nally , in Section 5 w e demonstrate an analysis of real daily currency data under b oth univ ariate and m ultiv ariate settings, and pro vide comparisons with existing metho ds. W e conclude with a discussion. 2 Lik eliho o d-fre e mo dels Computational pro cedures to sim ulate from p osterior distributions, π ( θ | y ) ∝ π ( y | θ ) π ( θ ), of parameters θ ∈ Θ giv en observ ed data y ∈ X , are w ell established (e.g. Bro oks et al. 201 0). Ho we ver when p oin twis e ev alua tion of the lik eliho o d f unction π ( y | θ ) is computa- tionally prohibitive or intractable, alternativ e pro cedures are required. Like liho o d- 4 free metho ds (also known as appr oximate B ayesian c omputation ) p ermit sim ulation from an appro ximate p osterior mo del while circum ve nting explicit ev aluation of the lik eliho o d function (T av a r ´ e et al. 199 7; Beaumon t et al. 2002; Marjo ram et al. 20 03; Sisson et al. 2 007; Rat mann et al. 2009). Assuming data simulation x ∼ π ( x | θ ) under the mo del giv en θ is e asily obtainable, lik eliho o d- free metho ds em b ed the p osterior π ( y | θ ) within an augmen ted mo del π LF ( θ , x | y ) ∝ π ǫ ( y | x, θ ) π ( x | θ ) π ( θ ) , (2.1) where x ∼ π ( x | θ ), x ∈ X , is an auxiliary para meter on the same space as the observ ed data y . The function π ǫ ( y | x, θ ) is t ypically a standard smo othing k ernel (e.g. Blum 20 0 9) with s cale parameter ǫ , whic h w eights the in tr a ctable posterior with high v alues in regions when the observ ed data y and auxiliary data x are simi- lar. F or example, unifo r m k ernels are commonplace in lik eliho o d-free mo dels (e.g. Marjoram et al. 20 0 3; Sisson et al. 20 07), although alternative s suc h as Epanec h- nik ov (Beaumont et al. 2002) and Gaussian kerne ls (Pe t ers et al. 2008) pro vide im- pro ve d efficiency . The resulting approx ima t io n to t he true p osterior targ et distribu- tion π LF ( θ | y ) ∝ Z X π ǫ ( y | x, θ ) π ( x | θ ) π ( θ ) dx = π ( θ ) E π ( x | θ ) [ π ǫ ( y | x, θ )] (2.2) impro ve s as ǫ decreases, and exactly recov ers the targ et p osterior as ǫ → 0, as then lim ǫ → 0 π ǫ ( y | x, θ ) b ecomes a p oin t mass at y = x (R eev es and P ettitt 2 005). P osterior sim ulation fro m π LF ( θ | y ) can then pro ceed via standar d simulation al- gorithms, replacing p oint wise ev aluations of π LF ( θ | y ) with Mon t e Carlo estimates through the expectation (2.2), based o n draws x 1 , . . . , x P ∼ π ( x | θ ) from the mo del (e.g. Marjora m et al. 2003). Alternatively , sim ulatio n from the join t p osterior π LF ( θ , x | y ) is a v ailable b y contriving to cancel the in tractable lik eliho o ds π ( x | θ ) in sample w eights or acceptance probabilities. F or example, imp orta nce sampling from the prior predic - 5 tiv e distribution π ( θ , x ) = π ( x | θ ) π ( θ ) results in an imp ortance w eigh t of π LF ( θ , x | y ) / π ( θ , x ) ∝ π ǫ ( y | x, θ ), whic h is free of lik eliho o d terms. See Sisson et al. (2009) f o r a discussion of margina l and j o in t-space lik eliho o d-free samplers. In general, the distribution of π ( x | θ ) will b e diffuse, unle ss x is discrete and dim ( x ) is small. Hence, generating x ∼ π ( · | θ ) with x ≈ y is improbable for realistic datasets y , and as a result the degree of computation required for a go o d lik eliho o d-free ap- pro ximation π LF ( θ | y ) ≈ π ( θ | y ) (i.e. with small ǫ ) will b e prohibitiv e. In practice, the function π ǫ ( y | x, θ ) is expressed through low dimensional v ectors of summary statis- tics, S ( · ), suc h tha t π ǫ ( y | x, θ ) weigh ts the in tractable p osterior through (2.1) with high v alues in regions where S ( y ) ≈ S ( x ). If S ( · ) is sufficie n t for θ , then letting ǫ → 0 r ecov ers π LF ( θ | y ) = π ( θ | y ) as b efore, but with more acceptable computationa l ov erheads, as dim( S ( x )) << dim( x ). As suf- ficien t summary statistics a r e generally unav ailable, the use of non-sufficien t statistics is commonplace. The effect of less efficien t estimators of θ in (2.2) is a mor e diffuse appro ximation of π ( θ | y ). Hence the c hoice of summary statistics in any application is critical, with the ideal b eing low-dimen sional, efficien t and near-sufficien t. In this article, w e implemen t the lik eliho o d-free sequen tial Mon te Carlo sampler of P eters et al. (2008), detailed in App endix A. As the class of particle-based algo- rithms is the most efficien t curren tly av ailable in lik eliho o d-free computation (e.g. McKinley et al. 200 9), and within this class, the sampler of P eters et al. (2008 ) is the only o ne to allow non-uniform functions π ǫ ( y | x, θ ), this sampler pro vides the b est com bination of efficien t simulation and flexible mo delling. 6 3 Ba y esian α -stable mo dels W e now dev elop univ ariate and m ultiv aria te Ba yes ia n α -stable mo dels. Unlik e exist- ing methods, lik eliho o d-free inf erence is indep enden t of mo del para meterization. 3.1 Univ ariate α -stable M o dels Denote the c hara cteristic function of n i.i.d. univ ar iate α - stable distributed random v ariables X 1 , . . . , X n b y Φ X ( t ) . A p opular and con ven ient parameterization is Φ X ( t ) =    exp  iδ t − γ α | t | α  1 + iβ tan π α 2 sgn ( t )  | γ t | 1 − α − 1  if α 6 = 1 exp  iδ t − γ | t |  1 + iβ 2 π sgn ( t ) ln ( γ | t | )  if α = 1 , (3.1) where sgn ( t ) = t | t | and i 2 = − 1 (e.g. Samoro dnitsky and T aqqu 1994). Many alter- nativ e parameterizations are detailed in Nola n (2007) and Z olotarev (198 6). Under (3.1), the in tractable stable densit y function is con tinuous and unimo dal, ta king sup- p ort o n ( −∞ , 0) if α < 1 , β = − 1; (0 , ∞ ) if α < 1 , β = 1 and ( −∞ , ∞ ) otherwise. Efficien t sim ulation of auxiliary data, x ∼ π ( x | θ ), under the mo del is critical for the p erformance of lik eliho o d-free metho ds (Section 2). Here, it is straigh tf orw ard to generate α - stable v ariates under the mo del defined b y the characteristic function (3.1) (e.g. Devroy e 1986; Nolan 200 7). This a pproac h is provided in App endix B. 3.1.1 Summary statistics A k ey comp o nen t of lik eliho o d-free inference is the av ailabilit y of lo w-dimensional, efficien t and near-sufficien t summary statistics. Since α -stable mo dels can p o ssess infinite v ariance ( α > 1) and infinite mean ( α < 1), this c hoice m ust b e made with care. Here w e presen t sev eral candidat e summary v ectors, S 1 – S 5 , previously utilized for parameter estimation in the univ ariate α -stable mo del. In Section 4 w e ev aluate 7 the p erformance of these v ectors, and pro vide informed recommendations for the c hoice of summary statistics under the lik eliho o d-f ree f ramew ork. S 1 McCul lo ch ’s Quantiles McCulloch (19 86) and McCullo c h (199 8) estimate mo del para meters based on sample quan tiles, while correcting for estimator sk ewness due to the ev aluation of b q p ( x ), the p th quan tile of x , with a finite sample. Here, the data x ( i ) are a rranged in ascending order a nd matc hed with b q s ( i ) ( x ), where s ( i ) = 2 i − 1 2 n . L inear interpolatio n to p from the tw o adjacen t s ( i ) v alues then establishes b q p ( x ) as a consisten t estimator of the true quantiles . In vers io n of t he functions b v α = b q 0 . 95 ( · ) − b q 0 . 05 ( · ) b q 0 . 75 ( · ) − b q 0 . 25 ( · ) , b v β = b q 0 . 95 ( · ) + b q 0 . 05 ( · ) − 2 b q 0 . 5 ( · ) b q 0 . 95 ( · ) − b q 0 . 05 ( · ) , b v γ = b q 0 . 75 ( · ) − b q 0 . 25 ( · ) γ then pro vides estimates of α, β and γ . Note that fro m a computational p ersp ectiv e, in v ersion of v α , v β or v γ is not required under like liho o d-fr ee methods. Finally , w e estimate δ b y the sample mean ¯ X (when α > 1). Hence S 1 ( x ) = ( b v α , b v β , b v γ , ¯ x ) . S 2 Zolotar ev ’ s T r ansforma tion Based on a transformatio n of data from the α -stable family X → Z , Zolotarev (1986) (p.16) pro vides an alternative para meterization of the α -stable mo del ( α, β , γ , δ ) ↔ ( ν, η , τ ) with a characteristic function of the form log Φ Z ( t ) = − exp n ν − 1 2 h log | t | + τ − i π 2 η sgn ( t ) i + C  ν − 1 2 − 1 o , (3.2) where C is Euler’s constan t, and where ν > 1 4 , | η | ≤ min { 1 , 2 √ ν − 1 } and | τ | < ∞ . This parameterization has the adv an ta g e that log arithmic moments ha v e simple expressions in terms of parameters to b e estimated. F or a fixed constant 0 < ξ ≤ 1 2 (Zolotarev 1 9 86 recommends ξ = 0 . 25) and fo r in t eger n/ 3, the tra nsfor mat ion is Z j = X 3 j − 2 − ξ X 3 j − 1 − (1 − ξ ) X 3 j , j = 1 , 2 , . . . , n/ 3 . 8 Defining V j = log | Z j | and U j = sgn( X j ), estimates for ν, η and τ are then given b y b ν = max { ˜ ν , (1 + | b η | ) 2 / 4 } , b η = E ( U ) , b τ = E ( V ) , where ˜ ν = 6 π 2 S 2 ( V ) − 3 2 S 2 ( U ) + 1, using sample v ariances S 2 ( V ) and S 2 ( U ). As b efore, δ is estimated by ¯ X (for α > 1), and so S 2 ( x ) = ( b ν , b η , b τ , ¯ x ). S 3 Pr ess’s Metho d O f Moments F or α 6 = 1 and unique ev aluation p oints t 1 , t 2 , t 3 , t 4 , the method of momen ts equations obtained f rom log Φ X ( t ) can b e solve d to obtain (Press 1972; W eron 2006) log( b γ ) = log | t 1 | lo g ( − log | Φ( t 2 ) | ) − log | t 2 | lo g ( − log | Φ( t 1 ) | ) log | t 1 /t 2 | , b α = log log | Φ( t 1 ) | log | Φ( t 2 ) | log | t 1 /t 2 | b β    b α, b γ = b u ( t 4 ) /t 4 − b u ( t 3 ) /t 3 ( | t 4 | b α − 1 − | t 3 | b α − 1 ) b γ b α tan  b απ 2  , b δ    b α = | t 4 | b α − 1 b u ( t 3 ) /t 3 − | t 3 | b α − 1 b u ( t 4 ) /t 4 | t 4 | b α − 1 − | t 3 | b α − 1 where b u ( t ) = tan − 1 [ P n i =1 cos( tx i ) / P n i =1 sin( tx i )]. W e a dopt t he ev a luation p oin ts t 1 = 0 . 2 , t 2 = 0 . 8 , t 3 = 0 . 1 and t 4 = 0 . 4 as recommende d b y Koutrouv elis (1980 ), and accordingly obtain S 3 ( x ) = ( b α, b β , b γ , b δ ) . S 4 Empiric al Char a cteristic F unction The empirical c haracteristic function, b Φ X ( t ) = 1 n P n j = 1 e itX j for t ∈ ( −∞ , ∞ ), can b e used as the basis for summary statistics when standard statistics are not a v ail- able. E.g. this may occur through the non- existenc e of moment generating functions. Hence, we sp ecify S 4 ( x ) = ( b Φ X ( t 1 ) , . . . , b Φ X ( t 20 )) where t i ∈ {± 0 . 5 , ± 1 , ± 1 . 5 , . . . , ± 5 } . S 5 Me an, Quantiles and Kolm o gor ov-Smirn ov Statistic The Kolmogorov-Smirno v statistic is defined as K S ( X ) = sup z | F X n ( z ) − F Y n ( z ) | , the largest absolute deviation b etw een the empirical cumulativ e distribution functions o f auxiliary ( X ) and observ ed ( Y ) data, where F X n ( z ) = 1 n P n i =1 I ( X i ≤ z ) and I ( X i ≤ z ) = 1 if X i ≤ z a nd 0 otherwise. W e sp ecify S 5 ( x ) = ( ¯ x, { b q p ( x ) } , K S ( x )), where the set of sample quan tiles { b q p ( x ) } is determined b y p ∈ { 0 . 01 , 0 . 05 , 0 . 1 , 0 . 1 5 , . . . , 0 . 9 , 0 . 95 , 0 . 99 } . 9 Lik eliho o d-free inference ma y b e implemen ted under an y parameterization whic h p ermits data generation under the mo del, and for whic h t he summary statistics a re w ell defined. F ro m the ab ov e S 1 – S 5 are join tly w ell defined for α > 1. Hence, to complete the sp ecification of the univ a riate α -stable mo del w e adopt the indep enden t uniform priors α ∼ U[1 . 1 , 2], β ∼ U[ − 1 , 1], γ ∼ U[0 , 300] and δ ∼ U[ − 300 , 300] (e.g. Buc kle 1995). Note that the prio r for α has a restricted domain, reflecting the use of sample momen ts in S 1 – S 3 and S 5 . F or S 4 w e ma y adopt α ∼ U(0 , 2]. 3.2 Multiv ariate α -stable Mo dels Ba ye sian mo del sp ecification and sim ulation in the m ultiv ariate α -stable setting is c hallenging (Nolan et al. 2001; Nolan 2 008; Samoro dnitsky and T a qqu 1 9 94). Here w e follow Nolan (200 8 ), who defines the m ultiv a r iate mo del f or the random v ector X = ( X 1 , . . . , X d ) ∈ R d through the functional equations σ α ( t ) = Z S d |h t , s i| α Γ ( d s ) ; β ( t ) = σ − α ( t ) Z S d sgn ( h t , s i ) |h t , s i| α Γ ( d s ) µ ( t ) =    h t , µ 0 i α 6 = 1 h t , µ 0 i − 2 π R S d h t , s i ln |h t , s i| Γ ( d s ) α = 1 where t = ( t 1 , . . . , t d ), s = ( s 1 , . . . , s d ), h t , s i = P d i =1 t i s i , S d denotes t he unit d - sphere, Γ ( d s ) denotes the unique spectral measure, and σ α ( t ) represen ts scale, β ( t ) sk ewness a nd µ ( t ) lo cation (through the v ector µ 0 = ( µ 0 1 , . . . , µ 0 d )). Scaling prop erties of the functional equations, e.g. µ ( r t ) = r µ ( t ), mean that it is sufficien t to consider them on the unit sphere (Nolan 1997). The cor r esp onding characteristic function is Φ X ( t ) = E exp ( i h X , t i ) = exp  − I X ( t ) + i  µ 0 , t  10 with I X ( t ) = R S d ψ α ( h t , s i ) Γ( d s ) , where the function ψ α is given b y ψ α ( u ) =    | u | α  1 − i sgn ( u ) tan  π α 2  α 6 = 1 | u | α  1 − i 2 π sgn ( u ) ln | u |  α = 1 . The sp ectral measure Γ( · ) and lo cation v ector µ 0 uniquely c haracterize the m ul- tiv ariate distribution (Samoro dnitsky and T aqqu 1994) and carry essen tial informa - tion relating to the dep endence b et we en the elemen ts of X . The contin uous sp ec- tral measure is typic a lly well approx imated b y a discrete set of k Dirac masses Γ ( · ) = k P j = 1 w j δ s j ( · ) (e.g. Byczk ow ski et al. 1993) where w j and δ s j ( · ) resp ectiv ely denote the weigh t and D irac mass of the j th sp ectral mass at lo cation s j ∈ S d . By simplifying the in tegral in I X ( t ), computation with the c haracteristic function Φ ∗ X ( t ) = exp {− P k j = 1 w j ψ α ( h t , s j i ) } b ecomes tra ctable and data g eneration from t he distribution defined by Φ ∗ X ( t ) is efficien t ( App endix B). As with the univ ar ia te case (3.1), standard parameterizations of Φ ∗ X ( t ) will b e discontin uous a t α = 1, result- ing in p o or estimates of lo cation and Γ( · ). In the multiv ariate setting this is ov er- come b y Z olotarev’s M-parameterization ( No lan et al. 2001; Nolan 2 0 08). Although lik eliho o d- free metho ds a re parameterization indep endent, it is sensible to w ork with mo dels with go o d lik eliho o d prop erties. In a Bay esian framew ork w e parameterize the mo del via the sp ectral mass, whic h in v o lv es estimation of t he w eigh ts w = ( w 1 , . . . , w k ) and locatio ns s 1: k ∈ S d × k of Γ( · ). F or k = 2 this correspo nds to s i = (cos φ i , sin φ i ) ∈ S 2 . More generally , for s i =  s 1 i , . . . , s d i  ∈ S d , w e use h yp erspherical co ordina t es φ i = ( φ i 1 , . . . , φ i d − 1 ), where φ i d − 1 = tan − 1  s d i s d − 1 i  , . . . , φ i 1 = tan − 1   q  s d i  2 +  s d − 1 i  2 + . . . + ( s 2 i ) 2 s 1 i   . W e define priors f o r the parameters ( w , φ 1: k , µ 0 , α ), with φ 1: k = ( φ 1 , . . . , φ k ), as w ∼ Diric hlet( s , . . . , s ), φ i j ′ ∼ U(0 , 2 π ), µ 0 j ∼ N( ξ , κ − 1 ), α ∼ U(0 , 2) fo r i = 1 , . . . , k , 11 j = 1 , . . . , d , j ′ = 1 , . . . , d − 1 and imp ose the ordering constrain t s 1 i − 1 ≤ s 1 i for all i o n the first elemen t of eac h v ector s i . Note tha t by treating the w eights a nd lo cations of the sp ectral masses as unkno wn para meters, these ma y b e iden tified with those regions of the sp ectral measure with significan t p o sterior mass. This differs with the approac h of Nolan (19 97) where the sp ectral mass is ev a luated at a lar g e num b er of deterministic grid lo cations. In estimating Γ( · ) | X a significan t reduction in t he n um b er of required pro jection v ectors is ac hieve d (Sections 3.2.1 a nd 4.2 ) . F urther, the ab o ve prio r sp ecification do es not p enalize placemen t of sp ectral masses in close pro ximity . While t his pro ve d adequate for the pr esen t ed analyses, alternative priors ma y usefully inhibit sp ectral masses a t similar lo cations (Piev atolo and Green 1998). 3.2.1 Summary statistics S 6 Nolan, Panorsk a & McCul lo c h Pr o je ction Metho d F or the d -v ariate α -stable observ at io ns, X i , i = 1 , . . . , n , w e tak e pro jections of X i on to a unit h yp ersphere in the direction t ∈ S d . This pro duces a set of n univ ariate v alues, X t = ( X t 1 , . . . , X t n ), where X t i = h X i , t i . The information in X t can then b e summarized b y an y of the univ ariate summary statistics S 1 ( X t ) , . . . , S 5 ( X t ). This pro cess is rep eated for m ultiple pro jections ov er t 1 , . . . , t τ . With the lo cation param- eter µ 0 estimated by ¯ x , fo r sufficien t n umbers of pro jection ve ctor s, τ , the summary statistics S 6 ( X ) = ( ¯ x , S s ( X t 1 ) , . . . , S s ( X t τ )) f or some s ∈ { 1 , 2 , 3 , 4 , 5 } , will capture m uc h of the informat io n con tained in the m ultiv aria te data, if S s is itself informative . The b est c hoice of univ ariate summary v ector S s will b e determined in Section 4.1. W e ado pt a randomized approac h to the selection of the pro jection v ectors t 1 , . . . , t τ , a voiding curse of dimensionalit y issues as the dimension d increases (Nolan 2 008). 12 4 Ev aluation o f mo del and sampler p er formance W e no w analyze the p erfo r ma nce of the Bay esian α -stable mo dels a nd lik eliho o d-free sampler in a sequence o f simulation studies. F or the univ ariate mo del, w e ev aluat e the capabilit y of the summary statistics S 1 , . . . , S 5 , and contrast the results with the samplers of Buckle (19 95) and Lombardi (2007). The p erformance of the m ultiv a riate mo del under the statistics S 6 is then considered for tw o and three dimensions. In the fo llo wing, w e implemen t the lik eliho o d-free seq uential Monte Carlo algo - rithm of P eters et al. (2008) (App endix A) in order to sim ula te from the lik eliho o d- free appro ximation to the true po sterior π LF ( θ | y ) ≈ π ( θ | y ) giv en by (2.2). This algorithm samples directly from π LF ( θ | y ) using densit y estimates based on P ≥ 1 Mon te Carlo draws x 1 , . . . , x P ∼ π ( x | θ ) from the mo del. W e define π ǫ ( y | x, θ ) as a Gaussian k ernel so that the summary statistics S ( y ) ∼ N ( S ( x ) , ǫ 2 Σ) for a suitably c hosen Σ. All inferences are based o n N = 1000 particles drawn from π LF ( θ | y ). Detail of algo rithm implemen tation is remov ed to App endix A for clarity of expo sition. 4.1 Univ ariate summary statistics and samplers W e sim ulate n = 200 observ atio ns, y , from a univ ariate α -stable distribution with parameter v alues α = 1 . 7, β = 0 . 9, γ = 10 and δ = 10 . W e then implemen t the lik eliho o d- free sampler targeting π LF ( θ | y ) for eac h of the univ aria t e summary statis- tics S 1 - S 5 described in Section 3.1.1, with uniform prior s for all parameters (Section 3.1). Alternativ e prior sp ecifications w ere inv estigated (Lomb ardi 2007; Nolan 1997), with little impact o n the results. P osterior minim um mean squared error (MMSE) estimates for eac h parameter, a ve r a ged ov er 10 sampler replicates are detailed in T able 1. Mon te Carlo standard errors are rep orted in parenthe ses. The results indicate that all summary v ectors apar t 13 from S 5 estimate α and δ parameters well, and for γ , S 3 and S 4 p erform p o orly . Only S 1 giv es reasonable results for β and for all parameters join tly . Figure 1 illustrates a progression of the MMSE estimates of eac h pa rameter using S 1 , from the lik eliho o d- free SMC sampler output for eac h sampler replicate. As the sampler progresses, the scale parameter ǫ decreases, and the MMSE estimates iden tify the true parameter v alues as the lik eliho o d-f r ee p osterior appro ximation improv es. The results in T able 1 are based on using P = 1 Mon te Carlo dra ws from the mo del to estimate π LF ( θ | y ) (c.f. 2.2) within the likelihoo d-free sampler. Rep eating the ab o ve study using P ∈ { 5 , 10 , 20 } pro duced v ery similar results, a nd so we adopt P = 1 for the sequel as the most computationally efficien t c hoice. F or comparison, w e also implemen t the auxiliary v ariable Gibbs sampler of Buc kle (1995 ) and the MCMC in ve rsion and series expansion sampler of Lombardi (2007), based on c hains of length 1 00,000 iterations (1 0 ,000 iterations burnin), and using their re- sp ectiv e prior sp ecifications. The Gibbs sampler p erformed p o orly for most parame- ters. The MCMC metho d p erformed b etter, but has larger standard errors than the lik eliho o d- free sampler using S 1 . The MCMC sampler (Lombardi 2007) p erforms lik eliho o d ev aluations via in ve r se F ast F ourier transform (F F T) with approxim ate tail ev aluation using Bergstrom ex- pansions. This appro ac h is sensitiv e to α , whic h determine s the threshold b etw een the FFT and the series expansion. F urt her, as the tail b ecomes fatter, a finer spacing of FFT abscissae is required to control the bias introduced outside of the Bergstrom series expansion, significan tly increasing computation. Ov erall, this sampler work ed reasonably for α close to 2, t ho ug h with deteriorating p erformance as α decreased. The G ibbs sampler (Buck le 19 9 5) p erformed extremely p o orly for most settings and datasets, ev en when using their prop osed change of v ariables transformations. As suc h, the results in T able 1 represen t sim ulations unde r whic h b oth Gibbs and MCMC 14 samplers p erformed credibly , thereby typifying their b est case scenario p erfo rmance. 4.2 Multiv ariate samplers W e consider v arying n umbers of discrete sp ectral masses, k , in the approximation to the sp ectral measure Γ ( · ) = k P j = 1 w j δ s j ( · ) . W e assume that the num b er of sp ec- tral masses is know n a priori , and denote the d -v ariate α -stable distribution b y S α ( d, k , w , φ 1: k , µ 0 ). Priors are sp ecified in Section 3.2. F ollo wing the analysis of Section 4.1, w e incorp orate S 1 within the summary vector S 6 . F or da t a sets of size n = 2 00, w e initially consider the p erformance of the biv ariate α -stable mo del, S α (2 , k , w , φ 1: k , 0 ), fo r k = 2 and 3 spectral masses, with resp ect to parameter estimation and the impact of the n um b er of pro jection vec to rs t 1 , . . . , t τ . The true and me a n MMSE estimates of each parameter, placing pro jection v ectors at the true lo cations of the sp ectral masse s, are presen ted in T able 2. In addition, results are detailed using τ = 2 , 5 , 10 and 2 0 randomly (uniformly) pla ced pro jection v ectors, in order to ev a luate the impact of sp ectral mass lo catio n uncertaint y . The lik eliho o d- free sampler output results in go o d MMSE parameter estimates, ev en for 2 r a ndomly placed pro jection vec tors. The parameter least accurately estimated is the lo cation v ector, µ 0 . Directly summarized b y a sample mean in S 6 ( · ), estimation of lo cation requires a la rge n umber of observ ations when t he data hav e hea vy tails. Figure 2 illustrates progressiv e sampler p erformance fo r the S α (2 , 2 , w , φ 1:2 , 0 ) mo del, with α = 1 . 7, w = ( π / 4 , π ) and φ 1:2 = ( π / 4 , π ). Eac h circular scatter plot presen ts MMSE estimates of w eigh t (radius) a nd angles (angle) of the t w o sp ectral masses, based o n 10 sampler replicates . The sequence of plots (a )–(d) illustrates the progression of the parameter estimates as the scale parameter ǫ (of π ǫ ( y | x, θ )) decreases (and hence the accuracy of the lik eliho o d- free a pproximation π LF ( θ | y ) ≈ π ( θ | y ) impro v es). As ǫ decreases , there is a clear mov emen t of the MMSE estimates 15 to wards the true angles and w eights, indicating appropriat e sampler p erformance. With sim ulated da t a sets of size n = 40 0 , w e extend the previous biv a riate study to 3 dimensions, with k = 2 discrete sp ectral masses. The true parameter v alues, and p osterior mean MMSE estimates and asso ciated standard errors, based on 10 sampler replicates, are presen ted in T able 3. Again, reasonable para meter estimates are obtained (giv en finite data), with lo cation ( µ 0 ) again the most difficult to estimate. In analogy with F igure 2 , progressiv e sample p erformance for the first sp ectral mass (with w 1 = 0 . 7 and φ 1 = ( π / 4 , π )) for decreasing scale parameter ǫ is sho wn in Figure 3. Based on 2 00 replicate MMSE estimates (for visualization purp oses), the shading o f each p oint indicates the v alue of w 1 as a p ercen t age (blac k=0%, white=100%), a nd the lo catio n on the sphere represen ts the angles φ 1 . F o r large ǫ , the MMSE estimates f or lo cation ar e unifo rmly distributed ov er the sphere, and the ass o ciated w eigh t ta k es the f ull range of p ossible v alues, 0 − 100%. As ǫ decreases, the es tima t es of spectral mass lo cation and w eight become strongly lo calized and cen- tered on the true parameter v a lues, a g ain indicating appropriate sampler p erformance. Similar ima g es are pro duced fo r the second discrete sp ectral mass. 5 Analysis o f exc hang e rate daily retu rns Our data consist of daily exc hange ra tes f or 5 differen t currencies recorded in GBP b et we en 1 Jan uary 20 0 5 and 1 D ecem b er 20 07. The data in volv e 1 065 daily-av eraged LIBOR (London interbank offered rate) observ ations y ′ 1 , . . . , y ′ 1065 . The standard lo g - transform generates a log returns series y t = ln  y ′ t +1 /y ′ t  . Cursory examination of eac h returns series rev eals clear non-Gaussian ta ils and/or sk ewness (T able 4, bo ttom). W e initially mo del eac h currency series as indep enden t dra ws from a univ ariate α - stable distribution. P osterior MMSE parameter estimates for eac h currency are giv en 16 in T able 4, based on 10 replicate lik eliho o d-free samplers using t he S 1 summary vec - tor. F o r comparison, w e also compute McCullo c h’s sample quan tile based estimates (deriv ed from S 1 , c.f. Sec t io n 3.1.1 ), and maxim um lik eliho o d estimates using J. P . Nolan’s ST ABLE progra m (av ailable online), using the direct searc h SPDF option with searc h domains g iv en by α ∈ (0 . 4 , 2], β ∈ [ − 1 , 1], γ ∈ [0 . 00001 , 1] and δ ∈ [ − 1 , 1]. Ov erall, t here is go o d agreemen t betw een Bay esian, lik eliho o d- and sample-based es- timators. All currency returns distributions are significan tly differen t from Gaussian ( α = 2 , β = 0), and exhibit similar fa mily parameter ( α ) estimates o v er this time p e- rio d. How eve r, t he GBP to YEN con v ersion demonstrates a significan tly asymmetry ( β ) compared to the o ther currencies. An in teresting difference b et w een the metho ds of estimation, is tha t McCulloch’s estimates of α differ considerably from the p osterior MMSE estimates, ev en tho ugh the latter are constructed using Mc Cullo c h’s estimates directly as s ummary statistics, S 1 . One reason that the Ba y esian estimates are more in line with the MLE’s, is that lik eliho o d- free metho ds largely ignore bias in estimators used as summary statistics (comparing the closeness b et ween biased o r un biased estimators will pro duce similar results – consider comparing sample and maxim um lik eliho o d estimators of v ar ia nce). The multiv ariate α -stable distribution assumes that its marginal distributions, whic h are also α -stable, p ossess iden tical shap e parameters. This prop ert y implies imp ortant pra ctical limitations, one of whic h is that it is only sens ible to join t ly mo del data with similar marginal shap e parameters. Accordingly , based on T able 4, w e now consider a biv ariate a nalysis of A UD a nd EUR O currencies . Restricting the analysis to the biv ariate setting a lso p ermits comparison with the biv ariate frequen tist a pproac h described in Nolan (1997) using the MVST ABLE soft ware (av ailable o nline). A summary illustration o f the discrete approx imatio ns to the underlying contin u- ous sp ectral mass is shown in Figure 4. Assuming k = 3 discrete spectral masses and 17 based on 10 lik eliho o d-free sampler replicates, the mean MMSE p osterior estimates (solid blac k line) with mean 3 σ p osterior credibilit y in terv a ls (dotted line), iden tify re- gions of high sp ectral mass lo cated at 2.7, 3.9 and 5 .6 , with resp ective we ig h ts 0.45 , 0 .2 and 0.35. Brok en lines in Figure 4 denote the frequen tist estimates of Nolan (1997), based on the iden tification of mass o ver an exhaustiv e mesh grid using 40 (dashed line) and 80 ( dash-dot line) presp ecified grid lo cations (pro jections). Ov erall, b oth appro a c hes pro duce comparable summary estimates o f the sp ectral mass appro ximation, although the lik eliho o d-free mo dels generate full p osterior distri- butions, compared to Nolan’s frequen tist estimates. The assumption of k = 3 disc rete sp ectral masses pro vides a parsimonious r epresen tatio n of the actual sp ectral mass. F or example, the sp ectral mass lo cated at 2 .7 accounts for the first tw o/three masses based on Nolan’s estimates (80/40 pro jections). While the frequen tist a ppro ac h is computationally r estricted to biv ariate inference, the lik eliho o d- free approa c h may naturally be applied in m uch higher dimensions. 6 Discuss ion Statistical inference for α - stable mo dels is challenging due t o the computational in- tractabilit y of the densit y function. In practice this limits t he rang e of mo dels fitted, to univ ariate and biv ariate cases. By adopting lik eliho o d-free Bay esian metho ds w e are able to circum v en t this difficult y , and pro vide approximate, but credible p osterior inference in the general multiv ariate case, at a mo derate computational cost. Critical to this approac h is the av ailabilit y o f infor mativ e summary statistics for the parame- ters. W e hav e show n that m ultiv ariate pro jections of data on to the unit h yp ersphere, in combination with sample quantile estimators, are adequate for this task. Ov erall, our approac h has a n umber of adv an tages o v er ex isting metho ds. There is 18 far greater sampler consistenc y than alternativ e samplers, suc h as the auxiliary Gibbs or MCMC in ve r sion plus series expansion samplers (Buc kle 1995; Lombardi 2007). It is largely indep enden t of the complexities of the v arious parameterizations o f t he α - stable c haracteristic function. The lik eliho o d- free approach is conceptually stra ig h t- forw ard, and scales simply and is easily implemen ted in higher dimensions (at a higher computational cost). Lastly , b y p ermitting a full Ba y esian m ultiv ariate analysis, the comp onen t lo catio ns and w eigh t s of a discrete approximation to the underlying con- tin uous sp ectral densit y are allow ed to iden tify those regions with highest p osterior densit y in a parsimonious manner. This is a considerable adv antage o v er highly com- putational frequen tist approac hes, whic h require explicit calculation of t he sp ectral mass ov er a deterministic and exhaustiv e g rid (e.g. Nolan 1997). Eac h analysis in this article used man y millions o f data-generations from the mo del. While computatio n fo r lik eliho o d-free metho ds increases with mo del dimen- sion and desired accuracy of the mo del appro ximation (through ǫ ), muc h o f this can b e o ffset through parallelization o f the lik eliho o d-f r ee sampler (P eters et a l. 20 08). Finally , while w e hav e largely fo cused on fitting α -stable mo dels in the lik eliho o d- free framew ork, extensions to mo del selection through Bay es factors or model a v erag - ing are immediate. One ob vious candidate in this setting is the unknow n n umber of discrete sp ectral masses, k , in the appro ximation to the con tinuous sp ectral densit y . Ac kno wledgmen ts YF and SAS are supp orted b y the Australia n Researc h Council Disco v ery Pro ject sc heme (DP08 77432 & DP109 2805). G WP is supp orted b y an AP A sc ho la rship and b y the School of Mathematics and Statistics, UNSW a nd CSIR O CMIS. W e thank M. Lom bardi for generously pro viding the use of his co de (Lombardi 2007), and M. Briers, S. G o dsill, Xiaolin Lou, P . Shev che nk o and R. W olp ert, for though tful discussions. 19 References Alder, R., R . F eldman, and M. S. T aqqu (1998 ). A pr actic al guide to he avy-tails: Statistic al te chniq ues fo r analysing he avy-taile d distributions . Birkh¨ auser. Beaumon t, M. A., W. Zha ng, and D. J. Balding (2002) . Approximate Ba yes ia n computation in p opulatio n genetics. Genetics 16 2 , 2025 – 2035. Bergstrom, H. (1953). O n some expansions of stable distributional functions. Ark. Mat. 2 , 375–37 8. Blum, M. G. B. (2009). Appro ximate Ba yes ian computation: a non- parametric p ersp ectiv e. T echnic a l rep ort, Univ ersit´ e Joseph F ourier, Grenoble, F rance. Bro oks, S. P ., A. Gelman, G. Jones, and X.-L. Meng (201 0 ). Handb o o k of Markov chain Monte Carlo . Chapman and Hall/CRC. Buc kle, D. J. (19 9 5). Bay esian inferenc e for stable distributions. Journal of the A m eric an Statistic al Asso ciation 90 , 605–613. Byczk owsk i, T., J. P . Nolan, and B. R a jput (1993). Approx imatio n of m ultidimen- sional stable densities. J. Multiv. Anal. 46 , 13–31. Casarin, R. (2004). Ba y esian inference for mixtures of stable distributions. Working p ap er No. 0428, CEREMADE, University Paris IX . Cham b ers, J., C. Mallows , a nd B. Stuc k (1976). A metho d for simulating stable random v ariables. J. Am. Stat. Asso c. 7 1 , 340–334. Correction (198 7), 8 2 , 704. Devro ye , L. (1986) . An automatic metho d f or generating r a ndom v ariat es with a giv en c haracteristic function. SIAM Journal o n Applie d Maths 46 , 698–7 19. Doganoglu, T. a nd S. Mittnik (199 8 ). An appro ximation pro cedure for asymmetric stable paretian densities. Computational Statistics 13 , 463–4 75. 20 DuMouc hel, W. (1975 ). Stable distributions in statistical inference: Information from stably distributed samples. J. A m. Stat. Asso c. 70 , 386–393. F ama , E. (19 6 5). The b ehaviour of sto c k mark et prices. J. Business 38 , 34–1 0 5. F ama , E. and R. Roll (1968 ) . Some prop erties of symmetric stable distributions. Journal of the Americ an Statistic al Asso cia tion 63 , 817 – 83. Go dsill, S. (1999 ). MCMC and EM-based metho ds for inference in hea vy-tailed pro- cesses with alpha stable innov ations. In Pr o c. IEEE Signal Pr o c e s s ing Works h op on Higher Or de r Statistics . Go dsill, S. (2000). Inference in symmetric alpha-stable noise using MCMC and the slice sample r . In Pr o c. I EEE I nternational Confer enc e on A c oustics, Sp e e ch and Signal Pr o c essing , V olume VI, pp. 3806–380 9 . Jiang, W. a nd B. T urn bull (2004). The indirect metho d: Inference based on in ter- mediate statistics – A syn thesis and examples. Statistic al Scienc e 19 , 238–2 6 3. Koutrouv elis, I. (198 0). R egr ession type estimation of the par ameters of stable laws . Journal of the Americ an Statistic al Asso cia tion 75 , 918 – 928. Kuruoglu, E. E., C. Molina, S. J. Go dsill, and W. J. Fitzgenrald ( 1997). A new ana- lytic represe n tation for the alpha stable pro babilit y densit y function. In H. Boz- dogan and R. So y er (Eds.), AMS Pr o c e e dings, Se ction o n Bayesian S tatistics . Levy , P . (1924). Theorie des erreurs. La loi de Gauss et les lois exce pt io nelles. Bul letin-So ciete Mathematique de F r anc e 52 , 49–85. Lom bardi, M. and S. Go dsill (20 06). On-line Bay esian estimation of AR signals in symmetric alpha-stable noise. IEEE T r a nsactions on Signal Pr o c ess ing . Lom bardi, M. J. (2007). Ba yes ian inference fo r alpha stable distributions: A r a n- dom walk MCMC approac h. Comp. Statist. and D ata Anal. 51 , 268 8 –2700. 21 Mandelbrot, B. (19 6 0). The P areto- L evey la w and the distribution of income. In- ternational Ec onomic R eview 1 , 79–106. Marjoram, P ., J. Molitor, V. Plag nol, a nd S. T av are (200 3). Mark ov c hain Monte Carlo without lik eliho o ds. Pr o c. Nat. A c ad. S ci. USA 100 , 15324–15 328. McCulloch, J. H. (1 986). Simple consisten t es t imators of stable distribution pa- rameters. Comm. Stat. Simulation and c omp utation. 15 , 1109– 1136. McCulloch, J. H. (1998 ). Numerical approximation of t he symme tric stable distri- bution and densit y . In E. Adler, R. R . F eldman, and T. M. Birkhauser (Eds.), A pr actic al guide to he avy tails : statistic al te chniq ues and applic ations. McKinley , T., A. R. Co ok, and R. Deardon (20 09). Inference in epidemic mo dels without likelihoo ds. The International Journal of Biostatistics 5: article 24. Melc hiori, M. R. (2006). T o ols for sampling multiv ariate archimedian copulas. www.YieldCurve.c om . Mittnik, S. and S. T. Rache v (1991 ) . Alernativ e m ultiv ariate stable distributions and their applications to financial mo delling. In E. Cambanis, G. Samoro dnit- sky , a nd M. S. T aqqu ( Eds.), Stable pr o c e s s es a n d r e la te d topics , pp. 107– 119. Birkhauser, Bo ston. Neal, R. (2003). Slice sampling. A nn als of Statistics 31 , 705 –767. Nikias, C. and M. Shao (1995). Signal pr o c es s ing with alpha stable distributions and appl i c ations . Wiley , New Y ork. Nolan, J. P . ( 1 997). Numerical computatio n of stable densities and distributions. Comm. Statisti. Sto chastic Mo dels 13 , 759 – 774. Nolan, J. P . (2007). Stable dis tributions: Mo dels for he avy-taile d da ta . Birkh¨ auser. 22 Nolan, J. P . (2 0 08). Stable distributions: Mo dels for hea vy-tailed data. T ec hnical rep ort, Math/Stat Department, American Univ ersit y . Nolan, J. P ., A. K. P anorsk a, and McCullo c h (2 001). Estimation of stable sp ec- tral measures, stable non-Gaussian mo dels in finance and econometrics. Math. Comput. Mo del ling 34 , 1113– 1122. P eters, G. W., Y. F an, and S. A. Sisson (2008). On sequen tial Mon te Carlo, partial rejection con trol and approximate Bay esian computation. T ec h. rep. UNSW. Piev atolo, A. and P . Green (1998). Boundary detection t hro ugh dynamic p o lygons. Journal of the R oyal Statistic al So ciety, B 60(3) , 609–626. Press, S. ( 1972). Estimation in univ ariate and m ultiv ariate stable distributions. Journal of the Americ al Statistic al Asso ciation 67 , 842 – 846. Ratmann, O., C. Andrieu, T. Hinkley , C. Wiuf, and S. Richards on ( 2 009). Mo del criticism based on lik eliho o d-free inference, with an application to protein net- w ork ev olution. Pr o c. Natl. A c ad. Sci. USA 106 , 1 0 576–1058 1 . Reev es, R. W. and A. N. Pettitt (2 0 05). A theoretical f r amew ork for appro ximate Ba ye sian computation. In A. R. F rancis, K. M. M a ta wie, A. Oshlac k, and G. K. Sm yth (Eds.), Pr o c. 20th In t. Works. S tat. Mo d., Austr alia, 2 0 05 , pp. 393 –396. Salas-Gonzalez, D., E. E. Kuruoglu, a nd D . P . Ruiz (200 6). Estimation of mixtures of symmetric alpha-stable distributions with an unknown n um b er of comp o- nen ts. IEEE I n t. Conf. on A c oustics, Sp e e ch and S i g . Pr o c. T oulouse, F r anc e . Samoro dnitsky , G. and M. S. T a qqu (1994). Stable non-Gaussian r andom pr o c esses: Sto chastic mo dels with infinite varian c e . Chapman and Hall/CR C. Sisson, S. A., Y. F an, and M. M. T anak a (2007). Sequen tial Mon t e Carlo without lik eliho o ds. Pr o c . Nat. A c ad . Sci. 104 , 1760– 1765. Errata (2009), 106, 16889. 23 Sisson, S. A., G. W. Peters , M. Briers, and Y. F an (20 09). Lik eliho o d-free samplers. T ec hnical rep ort, Unive r sity of New South W ales. T av ar ´ e, S., D . J. Balding, R. C. Griffiths, and P . Donnelly (1997) . Inferring coales- cence times from D NA sequence data. Genetics 145 , 505–518. W eron, R . (2006). Mo deling and for e c asting ele ctiricty lo ads and pric es: A statistic a l appr o ach . Wiley . Zolotarev, V. M. (198 6). One-Di m ensional Stabl e D istributions . T ranslations of Mathematical Monographs. American Mathematical So ciet y . 24 App endix A SMC sampler PRC-ABC algorithm (P et ers et al. 2008) Initialization: Set t = 1 and sp ecify tolerance sche dule ǫ 1 , . . . , ǫ T . F or i = 1 , . . . , N , sample θ ( i ) 1 ∼ π ( θ ), and set we ig h ts W 1 ( θ ( i ) 1 ) = π LF , 1 ( θ ( i ) 1 | y ) /π ( θ ( i ) 1 ). Resample: Resample N par ticles with resp ect to W t ( θ ( i ) t ) a nd set W t ( θ ( i ) t ) = 1 N , i = 1 , . . . , N . Mutation and correction: Set t = t + 1 and i = 1 : (a) Sample θ ( i ) t ∼ M t ( θ t ) and set w eight f o r θ ( i ) t to W t ( θ ( i ) t ) = π LF ,t ( θ ( i ) t | y ) / M t ( θ ( i ) t ). (b) With proba bilit y 1 − p ( i ) = 1 − min { 1 , W t ( θ ( i ) t ) /c t } , reject θ ( i ) t and go to (a ). (c) Otherw ise, a ccept θ ( i ) t and set W t ( θ ( i ) t ) = W t ( θ ( i ) t ) /p ( i ) . (d) Incremen t i = i + 1. If i ≤ N , go to (a). (e) If t < T then go to Resample. This algorithm samples N w eighted p articles from a sequence o f distributions π LF ,t ( θ | y ) giv en b y (2.2), where t indexes a sequence of scale parameters ǫ 1 ≥ . . . ≥ ǫ T . The final particles { ( W T ( θ ( i ) T ) , θ ( i ) T ) : i = 1 , . . . , N } , form a weigh ted sample from the targ et π LF ,T ( θ | y ) (e.g. P eters et al. 2008). The densities π LF ,t ( θ | y ) are esti- mated thro ug h the Monte Carlo estimate o f the expectation (2 .2 ) based on P draws x 1 , . . . , x P ∼ π ( x | θ ). F or the sim ulations presen ted w e implemen t the follow ing sp ecifications: for uni- v ariate α -stable mo dels θ = ( α, β , γ , δ ) and for m ultiv ariate mo dels θ = ( w , φ 1: k , µ 0 , α ); w e use N = 100 0 particles, initialized with samples from the prior; t he function π ǫ ( y | x, θ ) is defined b y S ( y ) ∼ N ( S ( x ) , ǫ 2 ˆ Σ) where ˆ Σ is an estimate of Cov( S ( x ) | ˆ θ ) based on 10 00 draws x 1 , . . . , x 1000 ∼ π ( x | ˆ θ ) giv en an appro ximate maxim um lik e- liho o d estimate ˆ θ of θ (Jia ng and T urn bull 2004); the mutation k ernel M t ( θ t ) = 25 P N i =1 W ( i ) t − 1 ( θ ( i ) t − 1 ) φ ( θ t ; θ ( i ) t − 1 , Λ) is a de nsity estimate of the previous particle po pula t io n { ( W t − 1 ( θ ( i ) t − 1 ) , θ ( i ) t − 1 ) : i = 1 , . . . , N } , with a Ga ussian k ernel density φ with cov ariance Λ; for univ ariate α - stable mo dels Λ = diag(0 . 25 , 0 . 25 , 1 , 1), and for m ultiv ariate mo d- els Λ = diag(1 , . . . , 1 , 0 . 25) (with D iric hlet prop osals a nd k ernel densit y substituted for w ); the sampler particle rejection threshold is adaptive ly determined a s the 90 th quan tile of the w eigh ts c t = b q 0 . 9 ( { W ( i ) t ( θ ( i ) t ) } ) , where { W ( i ) t ( θ ( i ) t ) } are t he N particle w eigh ts prior to particle rejection (steps (b) and (c)) at eac h sampler stage t (see P eters et al. 2008). F or eac h analysis w e implemen t 10 indep enden t samplers (in order t o monitor algo- rithm p erformance and Mon te Carlo v ariabilit y), each with the deterministic scale pa- rameter sequence: ǫ t ∈ { 1000 , 900 , . . . , 200 , 10 0 , 99 , . . . , 11 , 10 , 9 . 5 , 9 , . . . , 5 . 5 , 5 , 4 . 95 , . . . , 3 . 05 , 3 , 2 . 99 , 2 . 98 , . . . , 0 . 01 , 0 } . Ho w ev er, w e a da ptiv ely terminate all samplers at the largest ǫ v alue suc h that the effectiv e sample size (estimated b y [ P N i =1 [ W ( i ) t ( θ ( i ) t )] 2 ] − 1 ) consisten t ly drops b elo w 0 . 2 N ov er all replicate sampler implemen tations. App endix B: Data generation Simulation of univa ri a te α -stable data ( D uMouchel 1975, Chamb ers et al. 1976) 1. Sample W ∼ Ex p(1) to obtain w 2. Sample U ∼ Uniform[ − π / 2 , π / 2 ] to obtain u 3. Apply transformation to obtain sample y y =        S α,β sin α ( u + B α,β ) (cos u ) α/ 2  cos ( u − α ( u + B α,β )) w  1 − α α if α 6 = 1 2 π h  π 2 + β u  tan u − β ln π 2 w i cos u π 2 + β u i if α = 1 with S α,β =  1 + β 2 tan 2  π α 2  − 1 / 2 α and B α,β = 1 α arctan  β tan  π α 2  . In this case y will hav e distribution defined b y Φ X ( t ) with parameters ( α , β , 1 , 0). 26 4. Apply transformation to obtain sample y = γ y + δ with parameters ( α , β , γ , δ ). Simulation of d -dimensional, multivariate α -s tabl e data (Nolan 2007) 1. Generate Z 1 , ..., Z k i.i.d. random v ariables fro m the univ ar ia te α -stable distri- bution with parameters ( α, β , γ , δ ) = ( α , 1 , 1 , 0). 2. Apply the transformation Y =        k P j = 1 w 1 /α j Z j s j + µ 0 α 6 = 1 k P j = 1 w j  Z j + 2 π ln ( w j )  s j + µ 0 α = 1 with s 1 , . . . , s k , µ 0 ∈ S d . Note that while t he complexit y for generating realizations from a multiv ariate α -stable distribution is linear in the n umber of p oint masses ( k ) in the sp ectral represen t a tion per r ealizat io n, this metho d is strictly only exact for discrete sp ectral measures. Buc kle Lomb ar di S 1 S 2 S 3 S 4 S 5 α (1.7) 1.77 (0.18) 1.62 (0.10) 1.69 (0.06) 1.65 (0.07) 1.70 (0.06) 1.71 (0.04) 1.56 (0.05) β ( 0. 9) 0.54 (0.21) 0.86 (0.18) 0.86 (0.10) 0.65 (0.13) 0.31 (0.09) 0.38 (0.12) 0.49 (0.11) γ (1 0. 0) 18.17 (6.19) 9.59 (2.16) 9.79 (0.21) 10.44(0.56) 38.89 (6.34) 39.12 (5.92) 9.34 (0.14) δ (10.0) 12.30 (4.12) 9.70 (2.19) 10.64 (0.83) 9.31 (0.86) 10.25 (0.98) 10.83 (1.34) 11.18 (1.05) T able 1: Means and stand ard errors (in parentheses) of p osterior MMSE estimates of α , β , γ and δ under the univ ariate α -stable mo del, b ased on 10 sampler replicates. Paramete r v alues used for d ata simulatio n are given in the left column. Comparisons are b et ween the auxiliary v ariable Gibbs sampler metho d of Buc kle (199 5), the inv ersion MCMC metho d of Lom bard i (2007), and the lik eliho o d-free metho d, using summary statistic s S 1 – S 5 . 27 α µ 0 1 φ 1 φ 2 φ 3 w 1 w 2 w 3 T rue: k = 2 1.7 0 π 4 π – 0.6 0.4 – T rue: k = 3 1.7 0 π 4 π 3 π 2 0.3 0.25 0.45 k b α c µ 0 1 c φ 1 c φ 2 c φ 3 c w 1 c w 2 c w 3 Pro j ection vecto r s at locations of true spectral masses. 2 1.66 (0.04) 0.16 (0.19) 0.81 (0.65) 3.19 (0.46) – 0.55 (0.06) 0.45 (0.05) – 3 1.79 (0.02) 0.36 (0.18) 0.84 (0.27) 3.18 (0.29) 4.91 (0.24) 0.35 (0.05) 0.25 (0.04) 0.40 (0.05) 2 pro jection vect ors 2 1.67 (0.06) -0.13 (0.16) 0.73 (0.55) 3.58 (0.57) – 0.58 (0.09) 0.42 (0.10) – 3 1.76 (0.05) -0.16 (0.26) 0.91 (0.66) 3.65 (0.62) 4.85 (0.55) 0.36 (0.10) 0.24 (0.09) 0.40 (0.08) 5 pro jection vect ors 2 1.71 (0.05) 0.08 (0.14) 0.71 (0.60) 3.80 (0.67) – 0.60 (0.07) 0.40 (0.09) – 3 1.75 (0.04) 0.29 (0.17) 0.86 (0.62) 3.68 (0.52) 4.82 (0.41) 0.35 (0.09) 0.20 (0.07) 0.42 (0.09) 10 pro jection v ectors 2 1.72 (0.02) 0.21 (0.21) 0.75 (0.32) 3.31 (0.32) – 0.59 (0.05) 0.41 (0.07) – 3 1.73 (0.03) 0.25 (0.16) 0.76 (0.44) 3.31 (0.48) 4.78 (0.19) 0.34 (0.07) 0.24 (0.04) 0.42 (0.05) 20 pro jection v ectors 2 1.71 (0.03) -0.14 (0.14) 0.76 (0.36) 3.21 (0.23) – 0.63 (0.03) 0.37 (0.03) – 3 1.72 (0.02) 0.18 (0.23) 0.77 (0.32) 3.25 (0.31) 4.75 (0.15) 0.34 (0.04) 0.23 (0.03) 0.43 (0.03) T able 2 : Mean MMSE para meter estimates (and standard errors) fo r the biv ariate α -stable S α (2 , k , w , φ 1: k , µ 0 ) mo del, for k = 2 , 3 discrete sp ectral masses, calculated o ve r 10 replicate samplers. Pro jections v ectors are placed a t the true, and 2, 5, 10 and 20 randomly selected sp ectral mass lo cat io ns. The true v alue of µ 0 is the origin. α µ 0 1 φ 1 1 φ 1 2 φ 2 1 φ 2 2 w 1 w 2 T rue: k = 2 1.7 0 π 4 π π 2 3 π 2 0.3 0.7 k b α c µ 0 1 c φ 1 1 c φ 1 2 c φ 2 1 c φ 2 2 c w 1 c w 2 20 pro jection vec tors 2 1.71 (0.02) 0.53 (0.89) 1.12 (0.34) 3.81 (0.45) 1.84 (0.54) 4.24 (0.69) 0.28 (0.06) 0.72 (0.05) T able 3: Mean MMSE pa rameter estimates (and standard errors) for the triv ariate α -stable S α (3 , 2 , w , φ 1:2 , µ 0 ) mo del, with k = 2 discrete sp ectral masses, calculated o ve r 1 0 replicate samplers. The true v alue of µ 0 is the origin. 28 Currency Exc hange from GBP to AUD CNY EURO YEN USD b α 1.56 (0.03) 1.57 (0.02) 1.62 (0.04) 1.51 (0.04) 1.53 (0.02) Like l ihoo d b β 0.06 (0.03) 0.01 (0.009) -0.007 (0.08) -0.26 (0.09) -0.04 (0.03) free b γ 0.004 (4e-4) 0.003 (2e-4) 0.004 (1e-4) 0.003 (1e-4) 0.004 (3e-4) b δ 0.02 (0.01) 0.001 (0.0006) -0.03 (0.09) -0.06 (0.08) - 0.02 (0.07) b α 1.61 (0.05) 1.50 (0.05) 1.65 (0.05) 1.66 (0.04) 1.57 (0.05) MLE b β 0.08 (0.11) -0.01 (0.10) -0.10 (0.12) -0.46 (0.11) -0.01 (0.11) b γ 0.002 (7e-5) 0.002 (6e-5) 0.001 (4e-5) 0.002 (4e-5) 0.002 (1e-4) b δ -2e-4 (1e-4) -2e-5 (1e-4) 8e-5 (7e-5) 6e-4 (1e-4) 5e-5 (9e-5) McCullo c h’ s b α 1.39 1.38 1.47 1.38 1.39 quan tile b β 0.08 -0.003 -0.04 -0.18 0.001 estimates b γ 0.002 0.002 0.001 0. 002 0.002 b δ -4e-5 1e-6 1e-5 2e-4 5e-7 Kurtosis 8.39 9.11 15.60 6.29 4.98 Sk ewness 0.69 -0.42 -0.03 - 0.79 0.11 Std. dev. 0.004 0.004 0.003 0.004 0.003 Mean -4e-5 -4e-5 -9e-6 1e-4 7e-5 T able 4: P osterior MMSE estimates (Mon te Carlo errors) f rom the like liho o d-fr ee mo del, a nd maxim um lik eliho o d estimates (standard deviation). MLE’s, para meter estimates using McCullo ch’s quan tile (McCullo c h, 1998), and sample statistics (mean, standard deviation, sk ewness and kurtosis) obtained fr om J. P . Nolan’s ST ABLE soft ware, av ailable at academic2.am erican.edu/ ∼ jpnolan . 29 0 50 100 150 200 250 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 Likelihood−free sampler step (t) α 0 50 100 150 200 250 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Likelihood−free sampler step (t) β 0 50 100 150 200 250 0 10 20 30 40 50 60 Likelihood−free sampler step (t) γ 0 50 100 150 200 250 −15 −10 −5 0 5 10 15 20 Likelihood−free sampler step (t) δ Figure 1: T r aces of p osterior MMSE estimates of α , β , γ and δ under the un iv ariate α -stable mo del and lik eliho o d-fr ee sampler (summary statistics S 1 ), based on 10 sampler replicates. T races are sh o wn as a function ( x -axis) of sampler progression ( t ) and scale parameter redu ction ǫ t < ǫ t − 1 . Paramete r v alues used f or data generation are α = 1 . 7 , β = 0 . 9, γ = 10 and δ = 10. 30 0.5 1 30 210 60 240 90 270 120 300 150 330 180 0 (a) 0.5 1 30 210 60 240 90 270 120 300 150 330 180 0 (b) 0.5 1 30 210 60 240 90 270 120 300 150 330 180 0 (c) 0.5 1 30 210 60 240 90 270 120 300 150 330 180 0 (d) Figure 2: Circular scatter plot of MMSE estimates for k = 2 sp ectral mass angles (angle) and weigh ts ( radius) for biv ariate α - stable S α (2 , 2 , w , φ 1:2 , 0 ) mo del, with α = 1 . 7, w = (0 . 4 , 0 . 6) a nd φ 1:2 = ( π / 4 , π ). Plots (a)–(d) demonstrate ev olution of the estimates for decreasing scale parameter v alues ǫ , ba sed on 10 sampler replicates. 31 −1 −0.5 0 0.5 1 −1 0 1 −1 −0.5 0 0.5 1 ε = 100 −1 0 1 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 ε = 50 −1 0 1 −1 0 1 −1 −0.5 0 0.5 1 ε = 9 −1 0 1 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 ε = 4 −1 0 1 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 ε = 1.7 −1 0 1 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 ε = 1 0 10 20 30 40 50 60 70 80 90 100 Figure 3: Spherical heat map of MMSE estimates for the firs t of k = 2 discrete sp ectral masses, φ 1 , for the triv ariate α -stable S α (3 , 2 , w , φ 1:2 , 0 ) model. T rue v alues of the first sp ectral mass are w 1 = 0 . 7 (70%) and φ 1 = ( π / 4 , π ). P oint shading indicates MMSE v alue of w 1 as a p ercen tage. The plots demonstrate the ev olution of the estimates for decreasing scale parameter v alues ǫ , based on 200 sampler replicates. 32 0 1 2 3 4 5 6 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Location of spectral mass Normalised spectral mass Nolan (40 projections) Likelihood−free mean MMSE estimate Likelihood−free mean 3 σ C.I. Nolan (80 projections) Figure 4 : Es tima t es of s p ectral mass lo cation (x-axis) and cumulativ e w eight (y- axis) for AUD and EURO currencies data. Solid line denotes mean p osterior MMSE estimates of lik eliho o d-free SMC sampler output, and dotted line illustrates mean 3 σ p osterior credibilit y in terv a ls, based on k = 3 discrete spectral masses, 20 ran- domly placed pro jection v ectors and 10 replicate samplers. Broke n lines denote estimate of sp ectral mass using J. P . Nolan’s MVST ABLE softw are, a v ailable at academic2. american.ed u/ ∼ jpnolan , with (dashed line) 40 deterministic pro- jection lo cations and (dash-dot line) 80 pro jection lo cations. 33