A Bayesian Surrogate Model for Rapid Time Series Analysis and Application to Exoplanet Observations

Accepted for publication in Ba yesi an Analysis (2011) 6 , Number TBD, TBD A Ba y es ian Surrogate Mo del fo r Rapid Time Series Analysis and Application to Exoplanet Observations Eric B. F ord Univ ers ity of Florida Althea V. Mo orhead Univ ers ity of Florida Dimitri V eras Univ ers ity of Florida & Institute of Astrono m y Abstract. W e p resent a Bay esian surrogate mo del for the analysis of p eriod ic or quasi- p eriodic time s eries data. W e describe a compu t ationally efficient implementatio n that enables Bay esian mod el comparison. W e apply this mod el to sim ulated and real exoplanet observ ations. W e discuss th e results and demonstrate some of the chal lenges for applying our surrogate mo del to realistic exoplanet data sets. In particular, we find that analyses of real w orld data should pay car eful attention to the effects of uneven spacing of observ ations and th e choice of prior for the “jitter” parameter. Keyw ords: statistics: Ba yesian, mo del comparison, mo del selection, p erio dogram, numerical methods; exoplanets; observ ations: radial velocities, transit t iming v ari- ations 1 Motivation & Over view of Ex oplanet Observations Since the 1990’s , nearly 5 00 extras olar planets (or exoplanets) have b een discov ered around other sta rs in our g a laxy , yet o nly a few o f which have be e n observed dir ectly . In all other cas es, the planet’s pr esence ha s been inferre d fr o m its influence on the ligh t of its host star . 1.1 Doppler Observations The most pro ductive s uc h metho d to date has bee n observ ing the Do ppler shift of the star light due to the gravitational per turbations of the planet. F or a sing le planet on a cir cular orbit, the Doppler signatur e arises from a sinusoidal v ariation in the star’s velo c it y along the line of sig ht . More g e nerally , a single planet ca uses a p erio dic v ar ia tion in the s tellar velocity ( v ( t )) that follows the shap e predicted by Kepler’s laws of planetar y motion, which ca n usually be w ell-approximated b y the first few terms of a F ourier expansion (Konacki and Ma ciejewski 1999). If the amplitude of the fundamental term (with frequency f = 2 π /P , where P is the orbital per io d) in the F ourier expansion of v ( t ) is K 0 , then the co efficients of harmo nic ter ms (with frequency n f ) are o f order e n K 0 , where e is the orbital eccentricit y . Given the t ypical sig nal-to-nois e of detections c  2011 International So ciety for Bay esian Analysis http://ba.s tat.cmu.edu ba0001 2 Ba y esian Surrogate Mo del for Ti me Series Analysis and the typical exoplanet eccentricit y of less than 0.3 (and o ften less than 0.1 ), the use of a s few as t wo terms in this F ourier expansio n is often ac c urate to within observ ational uncertainties. If a star ha s multiple planets, its Doppler signature can b e muc h more complex. In some ca ses, the observed stellar velocities ( v k ) can b e well a ppr oximated as the line a r sup erp osition o f multiple planets on Keplerian orbits. In other cas es, the pla ne t- pla net interactions cause e ffects co mparable to the overall amplitude of the signal and thus must b e co nsidered. In a n y case, the observ ational sig nature is quasi- per io dic. A common question is whether a N p planet model is sufficient to describ e the av ailable observ ations o r whether the data demand at least N p + 1 planets. While an exploratio n of the full physical pa rameter space would b e computationa lly prohibitive, a low er dimensional surroga te mo del ca n b e quite useful for analyzing suc h a sys tem (V eras et al. 2 011). 1.2 T ransit Timing Va riations Recently , an extremely promising new metho d of detecting exo planets has bur st on to the scene. As a planet passes in front o f its host star (i.e., a tr ansit), the star’s brightness app ears to decrease. In an idealized tw o-b o dy system, the mid-times o f the transits ( t k ) are strictly p erio dic at the orbital perio d ( P ), so the k th trans it o ccurs at time t k = t 0 + k × P . If there are additiona l planets, then the times of the tra nsits deviate from a linear ephemeris (Agol et al. 2005; Holman and Mur ray 2005; F ord and Holman 2007). The p erturba tio ns of an additional planet can cause deviatio ns of the transit times that are simply sinusoidal in time, a p erio dic non-sinusoidal patter n or a very c o mplex q ua si- per io dic w av eform (Figur e 1, 2; F ord a nd Holman 2007; Nesvorn´ y and Mor bidelli 2008; Nesvorn´ y 2009; Nesvorn´ y and Beaug´ e 2010). While it is impractica l to explor e all the parameters o f a full ph ysica l mo del, a lower dimensional surro gate mo del may b e able to help identify reg ions of par ameter s pace that merit further inv estigation with a full ph ysica l mo del. 1.3 Relation to Previous Resea rch In this manuscript, we present a new method for ana lyzing p erio dic time ser ies data using a computationa lly efficient Bayesian surr ogate mo del. The deta ils of our mo del are chosen to facilita te the ana ly sis of exo planet observ ations. W e test our mo del by analyzing Doppler and tr ansit timing data sets. Thanks to the computational efficiency of o ur algo rithm, it is p ossible to apply it to a large libra r y of s imulated da ta se ts to understand how the mo del p erforms for different types o f planeta ry systems. Astronomers routinely apply Marko v chain Mo nt e Car lo (MCMC) techniques to per form Bay esian parameter estimation when analyzing Doppler obs erv a tio ns of an ex- oplanet host star (e.g., F or d 2005; Greg ory 2005). F or multiple planet sy s tems, MCMC metho ds a r e computationa lly intensive, e ven whe n the mo del ev a luation is p er formed ignoring the gravitational in teraction b e tw een the pla nets. While this is a go o d ap- proximation for a nalyzing the Doppler observ ations o f many systems, mut ual planetary int era c tions ca n be quite significant for planetary systems near a mean motion r esonance F o rd, Mo or head & V eras 3 (e.g., GJ 876; La ug hlin et al. 200 5). F ull n-b o dy integrations to account for thes e mu- tual planetar y interactions are muc h more computatio nally demanding. One appro ach is to develop and paralleliz e computatio nally efficient metho ds that allow one to use full n-b o dy integrations (e.g., Johnso n et al. 2011). While the GJ 876 data set has many high signal- to -noise o bserv a tions, most exoplanet host s ta rs hav e fewer o bs erv a tions and lo wer signal-to- no ise, resulting in weak constra int s o n many physical mo del pa ram- eters a nd making it even more difficult to sample these parameter s pa ces using MCMC. The us e of a lower-dimensional sur rogate mo del has the po tent ial to contribute to the analysis of such systems by identifying the perio dic ities that ar e statistically significant without introducing se veral a dditional parameters that a re po o rly constrained. Simi- lar metho ds a re routinely applied in a frequentist cont ext to identif y planet candidates from Doppler o bserv a tions (Cumming et al. 200 8). Our Bay esian surro gate mo del ca n be thought of as a Bayesian generaliz a tion o f the Lomb-Scargle p erio dog r am (Cumming 2004) that has b een further gener a lized to allow fo r m ultiple freque nc ie s , p erha ps due to the p ertur bations from additional pla ne ts or p erhaps due to significant eccen tric- it y (Kona cki and Maciejewski 1 999). Previously , a muc h more restricted version of the surrog ate mo del ( N f ,max = 1, N d,max = 0, see § 2.2) w as used to ev alua te strategies for scheduling Do ppler obs erv a tio ns (F ord 2008). The generalization in this manuscript allows for identifying multiple pe r io dicities, as is necessa ry for a pplication to eccentric and/or multiple planet sys tems. W e are optimistic that the surro gate mo del has ev en more p otential for con tributing to the analysis o f transit timing v aria tio n data. In the transit timing v ar iation method, the entire sig nal is due to the mutual gravitational p erturba tio ns. Given the highly non- linear na tur e o f the pr oblem (particular ly near r esonances), a ph ysica l mo del requir es per forming computationally ex pens ive n-b o dy in tegra tions. While it might b e practical to p erform MCMC sampling aro und o ne mode of the p osterior distribution while using full n-b o dy integrations, it is not practical to p erfor m a brute-force global sear ch of the high-dimensional parameter spac e while using full n-bo dy integrations (V e r as et al. 2011). The ev a luation of o ur surrogate model is orders of magnitude faster than an n-b o dy in tegra tion. Additionally , the surr ogate model is linear in most of its mo del parameters , allowing for efficien t identification of the mo des and integration over linear parameters , o nce we condition on the no n-linear par ameters. (W e p erform integration ov er non- linear parameter s via brute force, a s describ ed in the supplementary ma teri- als.) The sp eed of the surr ogate mo del makes it well-suited to exploring a br oad rang e of p os sible or bital co nfigurations. Once the surrog ate mo del ident ifies significa n t p eri- o dicities, n-bo dy in tegra tions can b e used to per form a more detailed exploration o f the full physical mo dels in r egions which could pro duce the p erio dicities identified by the surrog ate mo del. W e present results of o ur Bay esian surrog ate mo del applied to simu- lated transit timing data and discuss the implications of our res ults for the pr osp ects o f transit timing- based planet sea rches. 4 Ba y esian Surrogate Mo del for Ti me Series Analysis 2 Ba y esia n Surrog ate Mo del fo r Anal ysis o f quasi-P erio dic Time Ser ies 2.1 Data First, we describ e our general mo del for a na lysis of quasi- per io dic time serie s. In the case o f Doppler o bserv a tions, the independent v ariable, x would corresp ond to time ( t ) and y would cor resp ond to the star’s velocity ( v ( t )). In the case of transit timing observ ations, x would corres po nd to the transit num ber and y w ould cor resp ond to the mid- transit time. Each o bserv a tion ( y k ) is accompanied b y a n estimate of the measurement uncerta int y ( σ k ). The indep endent v ariables (i.e., transit num b er o r time of each o bserv a tion) are assumed to b e known precis ely . 2.2 Mo del W e ex plo re the use of a surrog ate model given by y ( x k ; θ ) = N f X i =1 [ S i sin (2 π f i x k ) + C i cos (2 π f i x k )] + N d X i =0 D i x i k , (1) where x k is the indep endent v a riable for the k th observ ation and and y ( x k ; θ ) is the predicted v alue of the observ able based on the surr ogate mo del with parameters θ . The surrog ate mo del parameters θ include: 1) the n umber of frequencies in the surrogate mo del ( N f ), 2 ) the frequencies ( f i ’s), 3) the amplitudes for those pe r io dicities ( S i ’s and C i ’s), 4) the order of the p olynomia l ter ms ( N d ), and 5) the p olynomial co efficients ( D i ’s). This mo de l exactly describ es time series which a re the sup erp osition o f p oly - nomial a nd sin usoidal signals. The surr ogate mo del can be used for B ayesian mo del compariso n to determine how complex a mo del (i.e., how many frequencies and/or po lynomial ter ms) is justified for a given data s et. In pr inciple, one could consider a lternative basis functions. W e fav or the use of sinusoids since they e xactly describ e the gravitational pe r turbation of a planet fo llowing a circular o rbit (and non-interacting with other planets). F urther , a sinusoid o ften provides a go o d firs t approximation to the per turbation of a planet on an eccentric or bit, given typical e c cent ricities and measurement precis ion. In practice, we will truncate the mo del to use just a few fr equencies and p olynomial terms, so as to pro vide an acceptable mo del for the observ a tions while facilita ting the efficient ev aluatio n of the mo del. If o bserv a tions s pan a s ufficiently long p erio d of time, then there is no need for p olynomial terms. In practice, the orbital p erio d of an outer planet (e.g., Saturn 29.5 years, Neptune 16 5 years) ca n b e muc h long er than the time span of obse r v ations (t ypically ∼ 1 - 10 years). If the or bital p erio d is muc h longer than the time span of observ ations than a simple linear term (constant r a dial acc e leration for Doppler data) is all that c an b e discerned from the av ailable data . In the case of a planet on a circular orbit, the “jerk” (i.e., deriv ative of acceler ation) b ecomes significant as the time span of observ ations approaches a qua rter of the orbital p erio d. F or the genera l case of a planet F o rd, Mo or head & V eras 5 on a n eccentric or bit, the jerk can b ecome eviden t after a g reater or lesser duration, depe nding up on the orientation of the orbit and phas e of the o bserv ations. In either ca s e, we cannot infer the physical par a meters with just tw o deriv atives measured, a s these ar e insufficient to ch ar a cterize a n orbit that r equres mor e pa rameters to fully sp ecify (sev en ph ysica l par ameters in genera l, and at least three even if we arbitr a rily as sume a circular, coplanar orbit). In these cases , using only one or tw o p olynomial terms dramatica lly increases the computational efficiency , since they elimina te the need to explor e a lar ge nu mber of fr e q uencies. Another a dv antage o f the p olynomial bas is is that a simple linear term corres po nds to the Doppler p e r turbation fr om a b o dy sufficiently distant that its orbital motion is insignificant over the time span of observ ations. A quadratic term corres p onds to physical motion of the p erturb er over the time span of observ ations. In ca ses where the time span o f observ ations is sufficient to discern a third p olynomial term, it b ecomes more difficult to justify a p olynomial mo del. On one hand, a cubic term is still muc h less computationally demanding than considering an additiona l harmo nic term (see § 2.5). On the other hand, acco rding to our physical mo del, a ll signals are per io dic on sufficiently long timescales and three deriv atives ar e sufficient to constra in the rang e o f plausible p erio ds and amplitudes. While three deriv atives are not sufficient to infer the shap e o f the wav eform, in many ca ses (e.g, a distant planet on nearly circular orbit) a simple s inusoid is a reas onable first approximation. Therefo re, we fav or considering a mo del with an additional ha rmonic term ov er including three o r more po lynomial ter ms. 2.3 Lik eliho o d W e take the likelihoo d of each observ a tion to b e L k ( x k , y k , σ k | θ ) = N ( y k − y ( x k ; θ ) , σ 2 k + σ 2 j ) (2) where N ( µ, σ 2 ) is a norma l dis tribution with mea n µ a nd v ariance σ 2 . Here θ =  N f , f 1 , ...f N f , S 1 , ...S N f , C 1 , ...C N f , N d , D 0 , ...D N d , σ j  . Note that we hav e expanded the se t of mo del pa rameters ( θ ) to include σ j , a “ jitter” parameter that is re la ted the amount o f scatter that is not acco unt ed for by the obs erv a tional uncertainties. The origin of the jitter need not be sp ecified. It could be due to inaccurate estima tes of the observ ational uncer taint ies or physical effects that a r e not included in a pa rticular mo del (e.g ., star sp ots, p -mo des , a dditional planets). W e a ssume that the obser v ational er rors a re uncor related in time, so the likelihoo d for a data set of N obs observ ations is given by L ( θ ) = N obs Y k =1 L k ( x k , y k , σ k | θ ) . (3) 6 Ba y esian Surrogate Mo del for Ti me Series Analysis 2.4 Prio rs In principle, o ne might aim to c ho ose priors that are based on the distributio n of mas s es and or bital p erio ds for exoplanets. Of course, characteriz ing those distributions is one of the primar y motiv ations to conduct ex oplanet s earches. While exoplanet searches hav e detected hundreds of planets , most of these have mass es compa r able to Sa turn or greater and orbital per io ds of a few years or less (Cumming et al. 2008). Inev itably , astronomer s pushing the frontiers of knowledge (e.g., searching for les s massive planets) will not know the intrinsic distr ibution of tho se pla nets’ orbital prop er ties . Thus, we choose broa d prior s ba sed on physical intuit ion and mathema tica l pr inc iples , as outlined in this se ction. F o r the num b er of frequencies in the sur roga te mo del, we adopt the following prior: p ( N f ) = (1 − 2 α + α N f ,max ) / (1 − α ) for N f = 0, p ( N f ) = α N f for 0 < N f ≤ N f ,max , and p ( N f ) = 0 for N f > N f ,max , where α parametrizes our prior b elief ab out the likelihoo d o f multiple freq uencies in the sig nal. The maximum num ber of frequencies to be c onsidered ( N f ,max ≥ 1) is chosen so as to provide eno ugh complexity to mo del the data while keeping the mo del ev aluation pr actical. The choice o f a geometric pro bability for an incr easing num b er of freq uencies has the adv antage that the prior ratio p ( N f = n + 1) / p ( N f = n ) is indep endent o f n (as long as n + 1 ≤ N f ,max ). If we a sso ciate each frequency with one planet and consider planet detection serially (i.e., first lo o k for evidence of one planet, next lo ok for evidence of another pla net), then the minimum Bay es factor to have a significant de tec tio n o f ea ch s uccessive planet is constant. The prior for the freq uencies themselves a re given by p (log f i ) ∼ U (log f min , lo g f max ), where U ( a, b ) is a unifor m distribution b etw ee n a and b . F or our applications, there are physical limits on the ra nge of viable freq ue ncie s. F or exa mple, f max could b e set by the shortest orbital p erio d in which a planet would b e a ble to survive for the age of the star , a nd f min could b e longe s t orbital per io d in which a pla net would b e able to remain b ound for the age of the star given p er tur bations fro m the ga lactic tidal field and passing star s. I n most ca ses, the time s pan of observ ations ( T obs ) will not b e s ufficien t to distinguish s uch long -p erio d sig nals fr om a low-order p olyno mia l. Th us, w e t ypically set f min ∼ 2 /T obs , as the sur rogate mo del is still able to mo del slow v ar iations and this aids in the rapid ev aluation of the surro gate mo del. Our choice of a pr ior that is flat in the lo garithm of the fr equency is motiv ated by the maximum entropy prior for a scale parameter . W e a lso ex per iment ed with a mo dified Jeffrey’s pr ior. W e found no significant difference in the results, a s long as there is a sig nificant detection. This can be under sto o d simply in terms of the c harac ter istic width of the peaks in the likeliho o d are muc h s maller than the domain of the frequencies. Thus, as long as the prior for frequency has significant supp ort across the entire doma in, the choice of the prior for frequency has minimal impact on the sha pe o f the p osterior (in lo cations with sig nifi- cant poster ior proba bilit y) unles s there is no t sig nificant evidence for any pe r io dicities. While we do no t attempt to justify o ur choice of prior based on the distributio n of exoplanet orbital p erio ds , we note that this a lter native approa ch would r esult in a fairly similar choice of prior , at least for freq ue ncie s b etw een ∼ 1/ (2 days) and ∼ 1/(20 00 days) (Cumming et a l. 200 8). F o rd, Mo or head & V eras 7 Each of the frequencies ( f i ) has an amplitude A i = p S 2 i + C 2 i and a phase φ i = atan2( − S i , C i ). W e choo se a uniform prior for each phase, p ( φ i ) = 1 / (2 π ). Physically , this co rresp onds to time inv ariance; i.e., the o ther planetary system do es not know what time we choos e to lab el a s t = 0. F or A i , the total a mplitude at frequency i , we a dopt a mo dified Jeffrey’s prio r , p ( A i ) ∼ (1 + A i / A o ) − 1 for 0 ≤ A i ≤ A max , where A max is the maximum plausible a mplitude. F or application to Doppler o bserv ations, A max could b e set by the maximum veloc ity p er tur bation by a planet (or binary star). F or application to trans it timing v ar iations, A max could b e set by the orbital p erio d. Note that, in pr inciple, A max could b e a function of f i . The other pa r ameter, A o , preven ts a divergence of the prior at sma ll amplitudes. F or so me applications A o may b e physically motiv ated. F or our applica tions, w e cho ose A o based o n the minimum detectable signal based on the av ailable da ta set, e.g., A o ∼  1 /σ 2 i  − 1 / 2 , where  1 /σ 2 i  − 1 / 2 is effectiv ely the weighted av era ge measurement precision. In practice, we find that any sufficiently small choice o f A o and large choice of A max gives similar re sults for pa rameter estimation for a given N f . The c hoice of A o , A max and α do affect the marginalize d pos terior probability ratio for N f . F o r our t ypical applications, the total signal amplitude is prop ortiona l to the mass of the plane t. While we do not attempt to justify our choice of prior bas e d on the distribution of ex oplanet ma s ses, we note that this a lter native appr oach would result in a fairly similar choice of prio r, at le a st for readily detectable planets (e.g., more massive than Saturn for Doppler sur veys). Present observ ations are only b eginning to provide significant constra in ts o n the distribution of planet ma sses for lo w-mass pla nets at s mall short orbital per io ds ( ≤ 50 days) (How ard et al. 20 10; Borucki et a l. 201 1; Y oudin 20 11). W e implement the ab ov e by using a “two-dimensional mo dified Jeffrey’s prio r ” for each pair of amplitudes S i and C i , p ( S i , C i ) = 1 /  4 π q S 2 i + C 2 i log (1 + A max / A o )  1 + q S 2 i + C 2 i / A o  (4) for q S 2 i + C 2 i ≤ A max , and p ( S i , C i ) = 0 for p S 2 i + C 2 i > A max . Inspection of Eqn. 4 shows that it is ess ent ially a mo dified Jeffery ’s prior for the to ta l a mplitude A i . F o r the num be r o f p olynomia l terms in the surrog ate mo del, we adopt a pr ior simila r to that for the num b e r of fr e q uencies: p ( N d ) = (1 − 2 β + β N d,max − 1 ) / (1 − β ) for N d = 1 , p ( N d ) = β N d − 1 for 1 < N d ≤ N d,max , and p ( N d ) = 0 for N d > N d,max , where β parametrizes our prior belief ab out the likelihoo d of higher order polyno mial ter ms being present in the signal. F or our applications, the maximum p olynomial order conside r ed ( N d,max ≥ 1) is typically s et to 1 and ra rely more than 2. The choice of a geometric probability for an increas ing num ber o f po lynomial terms has the adv a n tage that the prior ra tio p ( N d = n + 1) /p ( N d = n ) is indep endent o f n (as long as n + 1 ≤ N d,max ). F o r b oth the jitter and the polyno mial co efficients, i.e., each of B ∈ ( D i , σ j ), w e adopt a mo dified Jeffrey’s prior , p ( B ) = 1 / [2 | B | log (1 + B max /B o ) (1 + | B | / B o )] , | B | ≤ B max , (5) 8 Ba y esian Surrogate Mo del for Ti me Series Analysis where B o and B max are analog ous to A o and A max . W e adopt the same v alues as for A o and A max . F o r some data sets, it may b e wise to c ho ose a smaller B o , as the minim um measura ble magnitudes for the D i ’s depend o n the effective meas urement precision, the num b er o f o bserv a tions and the time-s pa n of obser v ations. In many cases, the amplitude of long-term trends is prop ortional to the p erturbing planet’s mass, motiv ating a prior p ( D i ) that is concentrated a t small signals a nd that has a similar s hap e to the prior for the amplitudes. W e a lso trie d using a unifor m prio r, p ( D i ) = 1 / (2 A max ) for | D i | ≤ A max . In practice, the likelihoo d is very sharply p eaked in D 1 , so results a re not sensitive to the choice of its pr ior. W e have only begun to exp eriment with the surr ogate mo del for N d = 2 (based on discussio n a t end of § 2 .2), and cautio n that further exp erimentation may b e needed for mo dels with N d ≥ 2. In this man uscript, we pr esent results based on a modified Jeffrey’s prior , as the jitter is a scale pa rameter and we hav e physical reason to limit p ( σ j ) at sma ll v alues of σ j , once it b ecomes sma ll co mpa red to the measurement pr ecision or astro physical effects tha t can cause non-gr avitational p er tur bations to the Doppler o r tra nsit timing signal (e.g., non- unifor m star s po ts). How ever, we found that in some cases our r esults could b e sensitive to the prior for the jitter (P ayne and F ord 2 011). F or e x ample, either impo sing a n upp er cutoff on the prior for σ j at ∼ 10m/ s or using a no r mal distribution for p ( σ j ) can significantly increase the poster ior probability for N f = n + 1 relative to that of N f = n or nar row the range of allow ed amplitudes. P revious studies o f the empirical distribution of Doppler jitter are bas ed o n r elatively small samples size s (W r ig ht 200 5), so they hav e little to s ay ab out the tails of the distribution. Based on our results, we encour age further obser v ations and statistical ana ly ses that could inform the choice o f prior for σ j . 2.5 Numerical Evaluation of Mo del and P osterio r PDF The surr ogate mode l was designed to provide a go o d a pproximation to Doppler or tran- sit timing observ ations and is lik ely to provide a reasona ble approximation of many other time se ries. A second importa nt feature of the sur roga te model is that it per - mits efficient e v aluatio n. Several tr icks to p erform an e fficien t brute force int egr ation are described in the supplemen tary material. A key feature of the mo del is that for given v alues of N f , N d , f i ’s and σ j , the mo del is linear. Thus, the integration ov er the remaining par ameters can b e pe rformed via linea r algebra and the La pla ce approxima- tion. Marginaliz ing ov er σ j and the f i ’s can b e p erformed using standar d numerical techn iques. By ev aluating the mo del co nditioned on N f and N d and margina lizing ov er the r emaining para meters, o ne can c o mpare the marginal p os terior pro babilities to ident ify v alues of N f and N d that pr ovide a go o d mo del for the observ ations without int ro ducing more mo del pa rameters than are justified given the av ailable data. F or our t ypical applica tions, the p os terior dominated by small v alues o f N f and N d , s o higher v alues need not b e cons ide r ed ex plicitly . Since the surrog ate mo del can b e in tegra ted nu- merically , it provides a quantitativ e bas is for Bay esian mo del compariso n and/or model selection. F o rd, Mo or head & V eras 9 3 Discussion of A pplication to Exoplanet Observ ations 3.1 Doppler Data One common challenge for exoplanet sea rches is deciding when the av ailable da ta pro- vides sufficient evidence to constitute a planet discovery . This is particula rly c hallenging in the case of stars with multiple planets, as the data can always b e b etter mo dele d by adding additional pla nets. The F ourier decomp osition of the Do ppler sig nature of a planet is dominated by p ower at the fr e q uency (1 /P ) c o rresp onding to the orbital pe - rio d ( P ). The p ower at the harmonic frequencies , 2 /P , 3 /P , ..., λ/P , dec r eases as e λ − 1 , where e is the orbital eccent ricity . The eccentricit y for a n ellipse is constrained to (0 , 1 ) and mo st known exoplanets have eccentricities smaller than 0 . 15. Thus, a pply ing the Bay esian surr ogate mo del to Doppler da ta sets is exp ected to result in a p os terior dis- tribution for the mos t significant frequency ( f 1 ) nea r one ov er the orbital p erio d of the planet which dominates the Doppler signature . The second most pro minen t freq uency ( f 2 ) could corres po nd to a harmonic of the first planet or to the o rbital per io d of a sec- ond planet. Some planetar y s y stems co nt ain tw o planets with o rbital p erio ds that differ by a factor of nea r ly 2 (e.g., Laughlin et al. 2005; W right et a l. 201 1; Lissauer et a l. 2011b). F or certain planet masses and eccentricities, the Doppler sig nature of t wo such planets can be mimick ed b y one planet with a mo re ec c en tric orbit (Giuppone et al. 2009; Anglada-E scud´ e et al. 201 0; Mo or head and F ord 2010). If there is actually only one planet, then one would ex pec t the next most pr ominent frequency to co r resp ond to twice the fundamen tal frequency ( f 2 = 2 × f 1 ) to within measurement precision. On the other hand, for systems that actually contain tw o planets, the sidere a l o rbital per io ds often deviate fro m exact resona nce (Lissauer et al. 2011 b). O ne might hop e that in cases wher e the ratio o f orbital per io ds differed from tw o, that the surrog ate mo del could reco gnize this difference ( δ ≡ f 2 − 2 × f 1 ), so one could infer that the observ ations were due to tw o plane ts , rather than one planet on a more highly eccen tric orbit (Giupp one et al. 20 09; Anglada- E scud´ e et al. 2 010; Mo or head and F ord 20 10). T o explo re this p ossibility , we applied the Bayesian surrog ate mo del to several s ets of Doppler obser v ations of exo planet host stars. W e fo cused on exo planets with a lar ge velocity amplitudes and b elieved to have a significa nt eccentricit y , as those systems provide the b est prosp ects for measuring the harmonic freq uencies prec isely . In partic- ular, we cho ose s ystems with K e 2 > 3 m/ s, where K is the velo city amplitude and e is the or bital eccentricit y . In all cases, the surroga te mo del efficiently found the funda- men tal frequency corres po nding to the orbital p erio d. As exp ected, f 1 is very tightly constrained and the ma r ginalized p os terior fo r additional fr equencies are significantly broader. In so me cases , we found that the marg inalized p osterior for f 2 did not co rre- sp ond to 2 × f 1 . In order to determine ho w often f 2 deviated from 2 × f 1 by chance, we p erformed a similar analysis on s everal simulated data sets. Each sim ulation was mo deled on simu- lated velocities that were ca lc ulated according to the b e s t-fit orbital per io d, amplitude, eccentricit y , arra ngement of p ericenter and o r bital phas e. W e a dded Gaussia n measure- men t noise with a scale set b y the cla imed measurement uncertaint y . The results for one 10 Ba y esian Surrogate Mo del for Ti me Series Analysis case (HD 162020 ) are shown in Figur es 3 & 4. If we use the same o bserv a tion times and uncertainties as the a ctual o bserv ations, then the marginal p osterio r distributions for P 2 = 1 /f 2 and P 1 / 2 = 2 /f 1 do not ov erlap. If we generate a s imilar data set, but with random obser v ation times, then the marginal p osterio r distributions for P 2 = 1 /f 2 and P 1 / 2 = 2 /f 1 do ov erla p. W e conclude that o ne must be very cautious o f alia sing due to unevenly spaced observ ations when interpreting the margina l p os terior distributions for f i ’s for r ealistic da ta sets. Our r esults sugg est that the problems of a lia sing co uld b e reduced by obtaining regular ly (or r andomly) spaced obse r v ations . Unfortunately , this is impractical for o b- serv a tions ma de from the surface o f the Ear th. Fir st, observ ations can not be made when the Sun is ab ove the hor izon (or even less than ∼ 12 ◦ below the hor izon) due to scattering of sunlight b y Earth’s atmospher e. As the Ea rth revolves around the Sun, the time of day at which a given sta r can b e obse rved from a given site changes. F or most stars, there are multiple months e ach year when high-precision Doppler o bserv ations are no t p ossible, since the star app e ars to o c lo se to the sun (after pro jecting onto the sky). Thus, for mos t stars , Doppler observ ations a re pr one to aliasing at frequencies as - so ciated with the sola r day and the sola r year. (In principle , obser v ations o f sta rs near the North or South p ole from an observ atory in the Arctic or An tarctic could av oid aliasing near the day . How ever, there are no obse r v ato r ies with high-pr ecision Doppler capability near either p ole due to lo gistical is s ues.) Second, while planet ho st stars are faint c ompared to the daytime sky brightness, they ar e m uch brighter than distant galaxies . Therefor e, time allo cation committees assign the v ast ma jor ity of obs erving time nea r new Mo on to extragala ctic as tronomers. Exoplanet searches a re typically a re assigned obser ving time near full Mo o n, intro ducing a liasing at fr e quencies asso ciated with the lunar month ( ≃ 29 . 5 days). In pra ctice, the b est wa y to avoid this is to dedicate an observ atory to Doppler o bserv a tions (e.g., the HARPS instrument a t the Europ ean Southern Obse rv ato ry’s 3.6 m telescop e in La Silla, Chile). Of course, this requires cons iderable r esources and is not an option for the world’s la rgest telescop es. Third, astr onomers often a ttempt to optimize the efficiency of their obs e rv atio ns o n a given night by obs erving each star when it is near the greatest (angular) a ltitude in the sky that night, a s this minimizes the amount of abso rption of starlight by the Earth’s atmosphere. This strategy intro duces yet another frequency asso ciated with the sider e al day (i.e., the time for a star to r e turn to the nearly same po int o n the sky from a given lo cation, roughly 23 hours, 56 min utes and 4 seconds). While these may seem like picky details, each of these timescales can be found in public data and in some cases c on- tribute to qua lita tive a m biguities in the or bital so lutio ns (Da wson a nd F abrycky 201 0). While atten tion to sch eduling can help impro ve efficiency of planet searches, ultimately weather (i.e., clo udy skies) will lead to data gaps and pr even t optimal exp erimental design for any Ear th-based o bserv a tory . Space- based observ atories are or ders of mag - nitude mor e ex pens ive to co nstruct and op era te. Therefore, we must develop to ols to analyze rea listic data sets. Studies such a s this will help us to use thos e to ols resp onsibly and reduce the r isk of making erroneo us claims. F o rd, Mo or head & V eras 11 3.2 T ransit Timing Va riations The r apidly incre asing num b er of known exoplanets that transit their host star ha s led se veral observers to measur e tr ansit times in hop es of detecting deviatio ns fr o m a linear ephemeris due to p er tur bations by another (p otentially non-transiting) planet (e.g., Miller -Ricci et al. 2008; Maciejewski et al. 2010). In most previous studies, the nu mber of precisely measured transit times has been insufficien t to fit a physical mo del. The pr osp ects for the tra ns it timing v ariation metho d are p o is ed to improve dra mat- ically thanks to NASA’s Kepler mission, whic h is obser ving over 10 0,000 sta rs near ly contin uously for 3 .5 years (Bo rucki et a l. 201 1; Steffen et al. 2010). Indeed, K epler re - cently rep orted the first strong evidence for transiting timing v ariatio ns in Ke ple r -9 and Kepler-11 sys tems (Holman et a l. 201 0; Lissauer et al. 2 0 11a). Interestingly , in b oth o f these cases, transit timing v ariatio ns were used to confirm planet candidates that had already b een iden tified via tr ansit and to constrain their masses and or bits. Kepler has als o identified dozens o f planet c andidates with putative transit timing v ariations (F o rd et al. 2011). Howev er, further o bserv a tions are needed in o r der to infer the mas s and or bital prop er ties of the p er tur bing bo dies. In principle, o ne could model the full light curve (i.e., observed brightness versus time; Carter et a l. 2 0 11). How ever, this would v a stly incr ease the computational time required and most of the information is contained in the tr ansit time. (Using data from NASA’s Kepler mission, there are t ypically ∼ 10 3 − 10 7 brightness measurements for each tr ansit time.) The transit time is the most precis ely measured par ameter for ea ch transit, and the transit times ar e s e ns itive to whether the transiting pla ne t is s lightly ahead or behind “ schedule” due to gravitational pe r turbations from other planets. The next b est-measured para meter is the transit dur ation whic h dep end on the orbit of the transiting pla net and the stellar ra dius (Mo orhea d et al. 201 1). The depth of each transit is primarily determined by the relative s izes of the planet and host star and is not affected by g ravitational p ertur bations o f o ther planets . The detailed shap e of e a ch transit also dep ends on detailed stellar prop erties (“limb darkening para meters”). F or a system of non-co planar planets, transit duration v ar iations may b e detectable. W e have fo cused our a nalysis in this pap er on coplana r systems viewed edge-on, s o as to reduce the dimensio nality of the para meter space to b e explored. While mo re parameter s are required to describ e no n-coplanar systems, in so me cases it may b e p ossible to der ive additional constraints on the orbits based on transit duration v ar ia tions, o r lack ther eof (Holman et al. 2010). Given the computatio na l cost of mo deling the full light curve, we recommend future r esearch to develop to o ls to analy ze a ser ies of transit times and transit dur ations. T o e xplore the potential for transit timing v aria tio ns to enable the detection of a non-transiting planet, we have genera ted ∼ 10 7 simulated data sets with a wide v a riety of orbital p erio ds, eccentricities a nd angles (V eras e t al. 2 011). W e apply the surrog ate mo del to s im ulated da ta sets to identify the dominant frequency a nd its amplitude. W e find that the surr ogate mo del can pr ovide an accurate mo del for some data s e ts, particularly those very near a mean motion reso nance, a regime which is par ticularly difficult to approximate ana lytically (Nesvorn´ y a nd Morbidelli 2008; Nesvorn´ y 20 0 9; 12 Ba y esian Surrogate Mo del for Ti me Series Analysis Nesvorn´ y and Bea ug´ e 2010). In other case s, the tr ansit timing v aria tion sig nature is more complex and w ould r equire sev era l frequencies to mo del adequately . F or extended data sets, this ca n be challenging, b oth due to computation time and a v a ilable memory . While the surro gate mo del can provide a reasona ble a pproximation to many tra n- sit timing signatures, the inferr ed mo del par ameters dep end sensitively on the orbita l phases, as well as more basic physical parameter s such as the planet mass a nd o r bital per io d. F urther, as one incr eases the num b er of obser v ations , the inferred par ameters often change significantly . This significa nt ly complicates the interpretation of the sur- rogate mode l outputs. While the inferre d mo del parameters will ev entually stabilize with a sufficient num b er and time s pa n o f o bserv ations, we find that infer red param- eters can contin ue to change e ven after s e veral years of obs erv a tions. W e c onclude that it may not alwa ys b e pra ctical to in vert transit timing v aria tions and infer the mass and orbita l pro per ties o f a non-tra nsiting pla net (Ragozzine a nd Holma n 2010). A mo re extended dis cussion of implications for transit timing planet sea rches is pre- sented separ ately (Pa yne et al. 2 0 10; V eras et al. 2 0 11). Of course , our r esults do not prov e that other ana ly sis techniques can not inv ert transit timing v ariatio ns. Ho wev er, our s ur roga te mo del was designed to ca ptur e the most imp ortant a s pec ts o f the prob- lem. Th us, our results are sug g estive that this problem may b e more general. Since the original submission of this pap er, the first confirmations of exo planets via the tra nsit timing v aria tion methods were published (Holman et al. 201 0; L is sauer et al. 2011 a). In these cases, e a ch of the detected planets transit the star, so transit timing v ar ia- tions were used to confirm planet candidates that had already b een identified by the standard transit technique. F ord et al. (2011) show ed that K e ple r can be exp ected to measure tra nsit timing v ariations for at least 12 systems with multiple transiting planet candidates. Based on a nalysis of the frequency of mult iple transiting pla net candidate systems (Lissauer et al. 201 1b), we exp ect that e ven more planets with transit timing v ar ia tions will b e s ig nificantly p erturb ed by a non- transiting pla net. Indeed, based o n the fir st four months of observ ations, Ke ple r has ident ified dozens of planet candidates with pr osp ective tra nsit timing v ariations, most in systems with only a sing le transiting planet candida te. Bo th our r esults and F ord et al. (2011) sugge s t that further observ a- tions will be neces sary b efore the mass es and or bits of putative additional pla nets can be deter mined. 3.3 Conclusions W e developed a Bayesian surrog ate mo del fo r analysis of time series data in g eneral and applied this mo del to tw o types o f exo pla net search data, Doppler and tra nsit timing v ar ia tions. The surro gate mo del can be ev aluated very r apidly and w e describ e a metho d for efficiently integrating over most mo del pa rameters. This a llows fo r calculating the (prop erly normalized) ma rginalized po s terior probability a nd assessing the p oster ior probability for a given p erio dicity . One strength of the surroga te mo del is for explor a tory data a nalysis. F o r exam- ple, as tr onomers ro utinely use the Lo mb-Scargle p erio do gram to sea rch Doppler data for a perio dic sig nature of a planet and to iden tify the range of perio ds that should F o rd, Mo or head & V eras 13 be explored with a more detailed mo del. O ne limiting case of the surr o gate mo del ( N f ,max = 1 and N d,max = 0) directly c orresp onds to the Bay esian gener alization o f the Lomb-Scargle p erio dogr am (Cumming 2 004). When a planet has a larg e eccentricit y o r one star hosts m ultiple pla ne ts , the Lo m b-Scar gle p erio dogr am typically reveals multi- ple sig nificant per io dicities. Previo usly , a stronomer s have dea lt with this by applying the Lo m b-Scar gle per io dogram to the residuals after subtracting the b est-fit sinusoidal or Kepler ian mo del. This approach ca n bias subsequent results, since the s ubtr acted mo del is not exact. F urther, assessing the significance of p eaks in the p erio dogr am of residuals is nontrivial. Mo st authors use a blind a pproach when sea rching for additional per io dicities, but others fav or using informatio n ab out the frequencies previously ident i- fied (e.g., Konacki and Maciejewsk i 1999; Dawson and F a brycky 2010). In pr a ctice, this can lead to a combersome dec ission tre e in a fr equentist con text. Our s urroga te mo del (with N f ,max > 1) represents a Bay e sian genera lization of iterative frequentist metho ds for a nalyzing perio do gram of residuals . The Ba yesian surrog ate mo de l provides a rigor - ous basis for calculating Bayes facto r of the mar ginalized p os terior probability for the nu mber of significa n t frequencies B n +1 ,n = p ( N f = n + 1 | x, y , σ ) /p ( N f = n | x, y , σ ). An- other adv a nt ag e of the s urroga te mo del is that by marg inalizing over the other mo del parameters (e.g., frequencies, a mplitudes, jitter, p olyno mia l ter ms), a spur ious false po sitive should b e less likely than a nalyzing r esiduals to only the b est-fit mo del. Finally , fo r many sys tems the surr ogate mo del ca n provide a low er-dimensional mo del that still captures the important (i.e., observ able) physical effects. F or e x ample, in a system of mult iple low-mass planets, a full physical mo del has a dimens io n of ≃ 7 N p , where N p is the n um b er of planets. If the system has planetary o r bits with small eccentricities and/or inclinations, then s everal of the mo del parameters ma y have no observ able effect. F or suc h systems, the surro gate model would b e able to descr ib e the sy stem a ccurately us ing a low er-dimensio nal parameter spa ce ( ≃ 3 N p ), g reatly increasing co mputational efficiency and p erha ps increa sing the sensitiv ity for detecting additional pla nets (due to the less extreme O ccam’s fa c to r). The surr ogate mo del is no t meant to replace o ther to ols for Bay esian par ameter estimation a nd mo del sele ction. It is still b eneficial to apply MCMC (and v a riants) for parameter e stimation using a mo r e physical mo del (e.g., F ord 2005; Greg o ry 2005; F ord 2006; Johnso n et a l. 201 1). Similarly , tools such as restricted Monte Carlo, imp or tance sampling and paralle l tempering can be helpful fo r calcula ting Bay es factor s using a more ph ysica l mo del (e.g ., F ord and Gre g ory 200 7; Gregor y 2 011). Since these to ols ar e m uch more computationally ex pens ive than the surr ogate mo del, they ar e most appropria te once a putative set of orbital perio ds has be en ide ntified (e.g ., by p erio do g ram analysis, surrog ate model, or human insp ection for sufficiently large s ignals). Thanks to the computatio nal efficiency of the surrogate mo del, we were able to ana- lyze numerous sim ulated da ta sets including multiple planet s y stems, something that would not hav e b een feas ible using previous ly av aila ble Bay esian metho ds such as MCMC. F or Doppler planet searches we find that realistic o bserving c adences ca n lead to significant alia s ing that prevents precisely testing whether there is a harmo nic rela tion- ship b etw een measured frequencies. Given the c lose relatio nship of our surro gate mo del to the widely used Lom b-Scar gle pe r io dogra m, our results also serve as caution regarding 14 Ba y esian Surrogate Mo del for Ti me Series Analysis the integration of results ba sed on p erio dogr am analyses (e.g., Konacki and Maciejewski 1999; Angla da-Escud´ e et al. 201 0; Da wson and F abrycky 2010). F o r analyzing transit timing v ariatio ns, we find the po sterior distribution for the surrog ate mo del para meters a re sens itive to the ex act orbital configura tio n. While the sensitivity to impo rtant ph ysical par ameters is adv antageous, sensitivity to parameter s that do not hav e dynamical significanc e makes in terpretation of the p os terior distri- bution for surroga te model parameter s more challenging. Unfortunately , the tra nsit timing sig nature often evolv es on a timescale co mpa rable to or longer than a r ealistic time spans for observ atio ns (e.g., 3.5-10 years for Kepler ). This makes it impractical to build a library o f p ossible tra nsit timing signatur es and the corre s po nding sur rogate mo del outputs. The surrogate mo del may still b e us eful for establishing the s ignificance of putativ e perio dic ities and/o r lo ng-term trends in transit timing data. In addition to the adv antages of a Bay esian approach, the computational efficiency o f the surro- gate mo del could b e useful for ana ly zing lar ge sets of s im ulated data s ets to aid in int erpr e tation of a tra nsit timing v ariation planet search. In both cas e s, we find that qualitative res ults (e.g., whether the Bay es fa ctor, B n,n +1 , for the sig nificance of a dditio nal frequency is greater or less than unity) can dep end o n the choice of pr ior for the jitter parameter ( σ j ). In this pap er we used a mathematically motiv ated prio r for σ j . Our result suggests that practical application of the surro gate mo del would significantly b enefit fr om further astronomical observ ations and statisti- cal a nalyses to asses s the empirical distribution of the stellar jitter for bo th Do ppler (e.g., W right 200 5) a nd transit timing observ atio ns (e.g., Holman et al. 201 0). Giv en the relationship of our sur roga te mo del to no n- Bay esian metho ds being employ ed by astronomer s, o ur results also sugg est ca ution in the in terpreta tion o f o ther results ba sed on frequentist a nalyses (which typically assume a single fixe d v alue o f the jitter). 4 Supplementa ry M aterial 4.1 Practical Mo del E valuation Integration over Linea r P arameters The surr ogate mode l is linear in the parameters: S i ’s, C i ’s, and D i ’s. F or a g iven choice of N f , N d , f i ’s a nd σ j , one c a n calculate the “ be st-fit” v a lues of S i ’s, C i ’s, and D i ’s via simple linear a lgebra. The sa me linear a lgebraic op erations allow the in tegr a ls ov er these linear para meters to b e e v aluated efficiently via the La place approximation (i.e., we e xpand the exp onent ab out the b est-fit para meters, keep the constant and second order terms, and extend the limits of in tegra tion to infinity; Cumming (2004)). F or our applications, the p oster ior probability is typically a smo o th function of σ j ’s, but can v ar y extremely rapidly with f i ’s. Th us, we recommend fix ing the choice of N f , N d , and f i ’s, a s that the int egr al ov er σ j can be ev aluated efficiently with standard numerical int egr a tion techniques. (In the sp e cial cas e of equal meas ur ement uncerta in ties, the int egr a ls can b e ca lc ula ted analytica lly .) F o rd, Mo or head & V eras 15 Brute F orce Integration over Frequencies Unfortunately , the in tegrals o ver the f i ’s m ust be ev aluated numerically . W e r ecom- mend discr etizing these in tegra ls and ev aluating them via brute fo r ce. The n umber of frequencies to b e ev aluated ( M ) can b e quite larg e (Cumming 2 004). F ortunately , there ar e several computational tricks that can sp eed up the calculation. In particular, most of the tr ig onometric functions can b e co mputed using trigono metric identities to improv e p erfor mance. W e find that brute force integration ov er one o r even tw o f i ’s is practical for realistic data sets. (F or larg e da ta sets, a la rge amount o f RAM may b e required for e fficie nt ev a luation.) Appro xima tions for Mo d els with Many Frequencies Unfortunately , brute force integration ov er more than tw o f i ’s ra pidly b ecome pro- hibitiv e. Thus, we intro duce the fo llowing approximation. W e start by p erforming a brute forc e ev aluation of the mo del conditioned on there b eing o ne frequency and then successively appr oximate the pos terior co nditioned on tw o frequencie s , three frequen- cies, etc.. When a pproximating the mo del conditioned on there b eing N f frequencies (and N f ≥ 2 ), we limit the s et of f i ’s with i < N f that are co nsidered to those which contributed significantly to the p osterior probability for the mo del conditioned on there being N f − 1 fr equencies. That is we set a thr eshold (e.g., ǫ = 1 0 − 4 ) and after ev aluating the mo del conditioned o n N f = 1 frequencies, we store the m 1 frequencies which hav e the larg est po sterior pr obabilities and collec tiv ely sum to a t least 1 − ǫ of the p os ter ior probability co nditioned on there b eing 1 frequency . At this p oint, we hav e searched M × (1 + m 1 ) frequencies (rather than M 2 ). Next, we estimate the marg inal poster ior probability for N f = 3 by co nsidering o nly tho se combinations of f 1 and f 2 that have the lar gest p osterio r probabilities and co llectively sum to at least 1 − ǫ o f the poste- rior proba bilit y conditio ne d on ther e b eing 2 frequencies. Thus, when ev aluating the mo del conditioned o n N f frequencies, there a r e no more than M Q N f − 1 i =1 m i frequencies to b e explicitly ev a luated, muc h less than M N f . In principle, this appr oximation can be r elaxed by pe rforming a full sea r ch over tw o freq ue nc ie s . In this case ev aluating the mo del conditioned o n N f frequencies, requir es no more than M 2 Q N f − 2 i =1 m i frequencies to b e explicitly ev aluated, which ma y or may not b e practical for a given data set. T run cating th e Num b er of Frequencies F o r data sets which a re well descr ib e d b e only one or a few frequencies, mo dels which include a large num b er of freq ue nc ie s will have a very small pos ter ior probability due to the Occam’s fa c tor as so ciated with the hig her-dimensional mo del. Thus, there is little po int in ev aluating mo dels with la rge N f . Th us, we recommend successively calc ula ting the pos ter ior pr obability co nditioned on N f ( p ( θ | N f , x k , y k , σ k )) and stopping at N f ,stop , such that p ( N f = N f ,stop | x k , y k , σ k ) ≪ P N f,stop − 1 i =0 p ( N f = i | x k , y k , σ k ) and appr oximate the r emaining mo dels a s p ( N f > N f ,stop | x k , y k , σ k ) ≈ 0. 16 Ba y esian Surrogate Mo del for Ti me Series Analysis Computat ional Cost The surrog ate mo del can b e ev alua ted muc h more quickly than an n-b o dy integration, explores a low er-dimensional pa r ameter space, and ta kes adv antage of the linea r ity of the mo del in most of the mo del parameters . Nevertheless, the str ong se ns itivit y to frequency dictates that w e m ust perform a fine sampling in frequency . F o r e x ample, consider a case of a ≃ 1 0 M E ar th -mass tra nsiting pla ne t with a n or bita l p erio d near 4 days and a s mall planet which is not observed to tra nsit with a p erio d near 8 days. With a mo dest eccentricit y (0.1), the transit timing v a riations o f the inner planet co uld be ∼ 10 minutes, compar able to the timing precision for each transit for a typical Kepler planet host star (F ord et al. 2011). Ov er 7.5 y ear s of observ ations (p ossible with an e xtended Ke ple r mission), o ne would observe roughly 680 trans its , allowing for an easy detec tio n o f such a single. If we ass ume a timing precision of 10 min utes and set f max to 1 days, then eac h in tegra l ov er fr equency requires consider ing M ∼ 18 0 , 000 frequencies. With a single core of a AMD O pteron 2 75 pr o cessor (2 .2 GHz), this takes ∼ 10, 16 or 37 seconds for models with N f = 1 , 2 or 3. Setting ǫ = 1 0 − 3 (see supplement ar y materia ls ), we needed to compute 1 ( N f = 1), 11 ( N f = 2) and 16 ( N f = 3) scans ov er frequency for b o th N d = 0 a nd 1, requiring a total of 12 C P U- min utes for one s ystem. Thus, the brute for c e ex ploration is relatively fast for a s ingle mo del. F or systems with no detectable signa l, the time required decrea ses significantly , as N f = 3 (or even 2) need not b e ex plored. On the other hand, the co mputation time p er system req uir ed g rows significantly as the signa l-to-noise increases , since the nu mber of freque ncie s sampled ( M ) mu st b e incr eased to avoid missing a narrow pea k in the p osterior density . O f cour s e, in these cases, the signa l is sufficiently la rge tha t fancy sta tistical metho ds are no t nece s sary to detect the dominant signal. The sp eed of the surr ogate mo del a llow ed us to analyze millions of simulated da ta sets and to e x plore the complex pa rameter space, using a cluster with hundreds o f AMD Opteron servers at the Universit y of Florida High-Performanc e Computing Center. References Agol, E., Steffen, J., Sari, R., a nd Clarkso n, W. (2005). “On detecting terres trial planets with timing o f gia nt planet trans its.” MNRAS, 35 9: 567– 579. Anglada-Es cud ´ e, G., L´ op ez-Mo rales, M., and Chamber s , J . E . (2 010). “How Eccentric Orbital Solutions Can Hide Planetar y Systems in 2:1 Resonant Orbits.” ApJ, 709: 168–1 78. Borucki, W. J ., Ko ch, D. G., Basri, G., Batalha, N., Brown, T. M., Bryson, S. T., Cald- well, D., Christensen-Dalsgaa rd, J ., Cochran, W. D., DeV o re, E., Dunham, E. W., Gautier, T. N., I I I, Geary, J. C., Gilliland, R., Gould, A., How ell, S. B., Jenkins, J. M., Latham, D. W., Lissauer, J. J., Mar cy, G. W., Row e , J., Sa sselov, D., B o ss, A., Charb onneau, D., Ciar di, D., Doyle, L., Dupree, A. K., F ord, E. B., F o rtney , J., Ho l- man, M. J ., Sea ger, S., Steffen, J. H., T arter, J., W elsh, W. F., Allen, C., Buchha ve, L. A., Christiansen, J. L., Cla rke, B. D., D´ es ert, J., E ndl, M., F abrycky , D., F ressin, F., Haas , M., Horch, E ., How ard, A., Isaacso n, H., K jeldsen, H., Ko lo dziejczak, J., F o rd, Mo or head & V eras 17 Kulesa, C., Li, J., Machalek, P ., McCarthy , D., MacQuee n, P ., Meib om, S., Miquel, T., Pr s a, A., Quinn, S. N., Quintana, E. V., Ragoz z ine, D., Sherr y, W., Shp or er, A., T enen baum, P ., T orres, G., Twicken, J. D., V a n Cleve, J ., and W alko wicz, L. (2011). “Charac ter istics of planetary candidates obser ved by Ke pler, I I: Analy sis of the first four months o f data.” ApJ, 73 6: 19–40. Carter, J. A., F a brycky, D. C., Ra gozzine, D., Ho lman, M. J., Quinn, S. N., Latha m, D. W., Buc hhav e, L. A., V an Cleve, J., Co chran, W. D., Co te, M. T., Endl, M., F ord, E. B., Haas, M. R., J enkins, J . M., Ko ch, D. G., Li, J., Lissa ue r , J. J ., Mac Queen, P . J ., Middour, C. K., O r osz, J. A., Rowe, J. F., Steffen, J. H., and W elsh, W. F. (2011). “KOI-12 6: A T riply E clipsing Hierarchical T riple with Tw o Low-Mass Stars.” Science , 3 31: 562 –. Cumming, A. (200 4). “Detectability of extraso lar planets in r adial velocity surveys.” MNRAS, 354: 1165 –1176 . Cumming, A., Butler, R. P ., Marcy, G. W., V ogt, S. S., W right, J . T., and Fischer, D. A. (2008 ). “The K e ck Planet Search: Detectability and the Minimum Mass and Orbital Period Distribution o f E x trasolar Planets.” P ASP, 120: 5 31–55 4. Dawson, R. I. and F abrycky, D. C. (2010 ). “Radial V elo cit y Planets De-alias ed: A New, Short Perio d for Sup er-E arth 55 Cnc e.” ApJ, 72 2: 937– 953. F o rd, E. B . (200 5). “Q uantif ying the Uncertaint y in the Orbits o f Extra solar Planets.” AJ, 129: 1706 –1717 . — (2 0 06). “Improving the Efficiency of Mar ko v Chain Monte Carlo for Analyzing the Orbits o f E xtrasola r P la nets.” ApJ, 642: 505 –522. — (200 8). “Adaptive Sc heduling Algorithms for Planet Sear ches.” AJ , 1 35: 10 08–10 20. F o rd, E . B. and Gre g ory, P . C. (20 07). “ Bay esian Mo del Selection and E xtrasola r Pla net Detection.” In G. J. Ba bu & E. D. F eigelson (ed.), Statistical Challenges in Mo dern Astronomy IV , volume 371 of Astrono mical Socie ty of the Pacific Conference Series , 189–+ . F o rd, E. B. and Holman, M. J. (2 007). “Using T r ansit Timing Observ ations to Sear ch for T ro jans of T ransiting Extra solar Planets.” ApJ, 664: L5 1–L54 . F o rd, E. B ., Row e, J. F., F abry cky, D. C., Ca rter, J., Holman, M. J., Lis s auer, J. J., Ragozzine, D., Steffen, J. H., B o rucki, W. J., Batalha, N. M., B r yson, S., Caldwell, D. A., Ga utier, T. N., I I I, Jenkins, J. M., K o ch, D. G., Li, J., Lucas, P ., Mar cy, G. W., Mullally, F. R., Quintana, E ., McCauliff, S., Thompson, S. E., T enenbaum, P ., Still, M., a nd Twicken, J. D. (201 1). “T ransit Timing Observ ations fro m K epler: I. Sta tistical Analysis of the First F o ur Months.” ApJS. Giuppo ne, C. A., T adeu dos Santos, M., Beaug´ e, C., F err az-Mello, S., and Mich tchenk o, T. A. (2009 ). “Detectability and E rror E stimation in Orbital Fits o f Resona n t E x - trasolar P la nets.” ApJ, 699: 132 1–133 2. 18 Ba y esian Surrogate Mo del for Ti me Series Analysis Gregor y , P . C. (20 05). “A Bay esian Analysis of Extrasola r Planet Data for HD 73 526.” ApJ, 6 31: 1198– 1214. — (2 011). “Bay esian Re-analys is of the Glies e 581 Exoplane t System.” MNRAS. Holman, M. J., F a brycky, D. C., Rago zzine, D., F or d, E. B., Steffen, J. H., W e ls h, W. F., Lissauer, J. J., Latham, D. W., Marcy, G. W., W a lko wicz, L. M., Batalha , N. M., Jenkins, J. M., Row e, J. F., Cochran, W. D., F ressin, F., T or res, G., Buchha ve, L. A., Sasselov, D. D., Borucki, W. J., K o ch, D. G., Basr i, G., Br own, T. M., Caldwell, D. A., Char bo nneau, D., Dunham, E. W., Ga utier, T. N., Geary , J. C., Gilliland, R. L., Haas, M. R., How ell, S. B., Ciardi, D. R., Endl, M., Fischer, D., F ¨ ur´ e sz, G., Hartman, J. D., Isaacson, H., Johnson, J. A., MacQ ueen, P . J., Mo orhe a d, A. V., Morehead, R. C., and Or osz, J. A. (20 10). “ K epler-9: A System of Multiple P lanets T ransiting a Sun-Like Star, Co nfirmed by Timing V ariations.” Science , 330: 5 1–. Holman, M. J. and Murray, N. W. (200 5). “The Use of T ransit Timing to Detect T errestrial-Ma ss Extrasola r Pla nets.” Science , 307: 128 8 –129 1 . How ard, A. W., Ma rcy, G. W., Johnson, J. A., Fischer, D. A., W right, J. T., Isaa cson, H., V alen ti, J. A., Anderson, J., Lin, D. N. C., a nd Ida , S. (2010 ). “The O ccurrence and Mass Distribution of Clo se-in Supe r -Ear ths, Neptunes, and Jupiters.” Science , 330: 65 3–. Johnson, J . A., Pa yne, M., How ard, A. W., Clubb, K. I., F or d, E. B., Bowler, B. P ., Henry, G. W., Fischer, D. A., Mar cy, G. W., Br ewer, J . M., Sc hw ab, C., Reffert, S., and Lowe, T. B. (20 11). “Retired a Sta rs and Their Companions. VI. A Pair of Int era cting Ex o planet P air s Around the Subgiants 24 Sextanis and HD 2 0096 4 .” AJ, 141: 16 –+. Konacki, M. and Maciejewski, A. J. (1999). “F requency Analysis of Reflex V elo cities of Stars with P la nets.” ApJ, 5 1 8: 442– 449. Laughlin, G., Butler, R. P ., Fis ch er, D. A., Marcy , G. W., V ogt, S. S., a nd W olf, A. S. (2005). “The GJ 876 Pla netary System: A Pr ogress Rep ort.” ApJ, 6 22: 1182– 1 190. Lissauer, J. J., F abrycky, D. C., F ord, E. B., Bo rucki, W. J ., F ress in, F., Marcy , G. W., Orosz, J. A., Rowe, J. F., T orr es, G., W elsh, W. F., Batalha, N. M., Br yson, S. T., Buchha ve, L . A., Caldwell, D. A., Car ter, J. A., Charbo nneau, D., Christiansen, J. L., Co chran, W. D., Desert, J., Dunham, E. W., F anelli, M. N., F ortney, J. J., Gautier, T. N., I I I, Gear y, J . C., Gillila nd, R. L., Haas, M. R., Hall, J. R., Holman, M. J., K o ch, D. G., Latham, D. W., Lo pez , E ., McCauliff, S., Miller , N., Mor ehead, R. C., Q uintana, E . V., Rag o zzine, D., Sa sselov, D., Short, D. R., and Steffen, J. H. (2011a ). “A closely pack ed system of low-mass, low-densit y planets transiting Keple r - 11.” Nature, 470 : 53–58. Lissauer, J. J., Ragozzine, D., F abrycky, D. C., Steffen, J. H., F o rd, E. B., Jenkins, J. M., Shpo rer, A., Holma n, M. J., Row e, J. F., Quintana, E . V., Ba talha, N. M., Bo rucki, W. J., Bryso n, S. T., Caldwell, D. A., Carter, J. A., C ia rdi, D., Dunham, E. W., F o rd, Mo or head & V eras 19 F o rtney, J. J., Gautier, T. N., I I I, Howell, S., Ko ch, D. G., Latham, D. W., Mar cy, G. W., Mo rehead, R. C., a nd Sasselov, D. (20 11b). “Architecture and Dynamics of Kepler’s Candidate Multiple T ransiting Planet Systems.” ArXiv e-pr int s . Maciejewski, G., Dimitrov, D., Neuhaeus e r , R., Niedzie ls ki, A., Rae tz, S., Ginski, C., Adam, C., Ma rk a, C., Moualla, M., and Mugrauer , M. (2 0 10). “T ransit timing v ar i- ation in exo planet W ASP-3b.” MNRAS, 407 : 2625– 2 631. Miller-Ricci, E., Row e, J. F., Sasselov, D., Matthews, J. M., Guenther, D. B., Kuschnig, R., Moffat, A. F. J., Rucinski, S. M., W alker, G. A. H., and W eiss, W. W. (2008 ). “MOST Space-ba s ed Pho tometry of the T ransiting Exoplanet System HD 209 458: T ransit Timing to Search for Additiona l Planets.” ApJ, 6 82: 586–5 9 2. Mo orhead, A. and F ord, E . (2010 ). “Resolving the degenerac y b etw een eccentric planets and 2:1 mean motion reso nances.” In K. Go ˙ zdziewski, A. Niedzielski, & J. Schneider (ed.), E AS P ublica tions Ser ies , volume 42 of EAS Publications Series , 16 1–164 . Mo orhead, A. V., F ord, E. B., Morehead, R. C., Row e, J., Borucki, W. J ., Ba talha, N. M., Bryson, S. T., Caldw ell, D. A., F abr ycky, D. C., Gautier, T. N., II I, Koch, D. G., Holma n, M. J., Jenkins, J. M., Li, J ., Lissa uer, J. J., Lucas, P ., Marcy , G. W., Quinn, S. N., Quintana, E., Ragozzine, D., Shpor er, A., Still, M., and T orr es, G. (2011). “The Distribution of T ransit Durations for Kepler Planet Candidates and Implications for their O rbital Ec cent ricities.” ApJS. Nesvorn´ y, D. (2009). “T ransit Timing V ariatio ns for Eccentric and Inclined Exoplanets.” ApJ, 7 01: 1116– 1122. Nesvorn´ y, D. a nd Bea ug´ e, C. (2010 ). “F ast Inversion Method for Determination of Planetary Parameters from T ransit Timing V a riations.” ApJ, 709 : L44–L48. Nesvorn´ y, D. a nd Mor bidelli, A. (2008 ). “ Ma ss a nd O rbit Deter mination from T ransit Timing V a riations of Ex oplanets.” ApJ, 68 8: 636– 646. Pa yne, M. J. and F ord, E. B. (2011). “An Analysis o f J itter and T r ansit Timing V a riations in the HA T-P-13 System.” ApJ, 729: 98–+ . Pa yne, M. J., F or d, E. B., and V eras, D. (2010 ). “T ransit Timing V ar iations for Inclined and Retro g rade Exo planetary Systems.” ApJ, 712 : L86–L92 . Ragozzine, D. a nd Holman, M. J. (20 10). “The V alue o f Systems with Multiple T ran- siting Pla nets.” ArXiv e-prints . Slawson, R. W., Pr sa, A., W elsh, W. F., Oros z, J. A., Ruck er, M., Ba talha, N. M., Doyle, L. R., E ngle, S. G., Co nroy, K., Coughlin, J ., Ames Gregg , T., F etherolf, T., Short, D. R., Windmiller, G., F abrycky, D. C., How ell, S. B ., J enkins, J. M., Uddin, K., Mullally, F., Seader, S. E., Thompson, S. E., Sanderfer, D. T., Borucki, W., and Ko ch, D. (2011 ). “K e ple r Eclips ing Binary Stars. I I. 2165 Eclipsing Binaries in the Second Data Relea se.” ArXiv e - prints . 20 Ba y esian Surrogate Mo del for Ti me Series Analysis Steffen, J. H. and Agol, E. (200 5). “ An analysis of the tra ns it times of T rES-1b.” MNRAS, 364: L96–L1 00. Steffen, J. H., Bata lha, N. M., Bo rucki, W. J., Buc hhav e, L. A., Caldwell, D. A., Co chran, W. D., Endl, M., F abrycky, D. C., F ressin, F., F ord, E. B ., F ortney, J. J., Haas, M. J ., Holma n, M. J., Is aacson, H., Jenkins, J . M., K o ch, D., La tha m, D. W., Lissauer, J. J., Mo orhead, A. V., Morehe a d, R. C., Mar cy, G., Mac Q ueen, P . J., Quinn, S. N., Ra g ozzine, D., Rowe, J. F., Sasse lov, D. D., Seage r , S., T orres, G., and W elsh, W. F. (2010). “ Five Kepler target stars that show multiple transiting exoplanet candidates.” ApJ, 725: 12 26–1 2 41. Steffen, J. H., Quinn, S. N., Borucki, W. J., Bruga my e r, E., Br y son, S. T ., Buchha ve, L. A., Co chran, W. D., Endl, M., F a br ycky , D. C., F or d, E. B., Ho lma n, M. J., Jenk- ins, J., K o ch, D., Latha m, D. W., Ma cQueen, P ., Mullally , F., Prs a, A., Ragozz ine , D., Row e, J . F., Sander fer, D. T., Seader, S. E., Short, D., Shp orer, A., Tho mpson, S. E., T orr es, G., Twick en, J. D., W els h, W. F., and Windmiller, G. (2011). “ The architecture of the hiera rchical triple s tar KOI 928 from eclipse timing v ariatio ns seen in K epler photo metry .” ArXiv e-prints . V er as, D., F ord, E . B., and Pa yne, M. J. (2 0 11). “Quantifying the Challenges of De- tecting Unseen Pla neta ry Companions with T ransit Timing V ariations.” ApJ, 72 7: 74–+. W r ight, J. T. (2005 ). “Radia l V elo city Jitter in Star s from the Califo rnia and Car negie Planet Sea rch at Keck Obser v ator y.” P ASP, 117: 657 – 664. W r ight, J . T., V er as, D., F ord, E. B., Johns on, J. A., Marc y , G. W., How ard, A. W., Isaacso n, H., Fischer, D. A., Spro nck, J., Anderson, J., and V alenti, J. (201 1). “ The California Planet Surv ey . I II. A Possible 2 :1 Resonance in the Exoplanetar y T riple System HD 37 124.” ApJ, 7 30: 93–+. Y o udin, A. N. (2011 ). “The Exo planet Ce nsus: A General Metho d, Applied to Kepler .” ArXiv e-prints . Ab out the Authors Eric B. F ord is an associate professor of astronom y at the Universit y of Flori da. His re- searc h fo cuses on extrasolar p lanets, includ ing b oth theory and the statistical interpretation of exoplanet observ ations. A lthea V. Mo orhead is a p ostdo ctoral associate study ing planet formation and contributing to NAS A’s K epler mission. Dr. Mo orhead applied th e surrogate mod el t o real and simulated radial velocity observ ations. Dimitri V eras w as a p ostdo ctoral associate stu dying planetary orbital dyn amics and the transit timing meth od. Dr. V eras ap- plied the surrogate mo del to sim ulated transit timing v ariations data sets. Dr. V eras is no w a p ostdoctoral asso ciate at the Institute of Astronomy in Cambridge, UK. Acknowledgments F o rd, Mo or head & V eras 21 The auth ors wish to th an k Mala y Ghosh, Phil Gregory , T om Loredo and Matthew Pa yne for v aluable d iscussions and feedback. This material is based up on work supp orted by the National Science F oundation under Grant No. 0707203. Additional supp ort for this work was p ro vided by th e National Aeronautics and Space Administration under gran t NNX09AB35Gb issued through the Origins of Solar Systems p rogram and grant NNX08AR04G issued through the Kepler P articipating Scientist Program. 22 Ba y esian Surrogate Mo del for Ti me Series Analysis Figure 1 : In this figure, we show the tr ansit timing signature for a ser ies of tw o pla net systems, ea ch extending for 3 .5 years, the nominal mission lifetime for NASA’s K epler space obser v ator y . In ea ch case, transit times are for a 20 Ear th-mass planet following an initially cir cular o rbit with a star -planet separa tion of approximately 0.05AU. (One A U is the av erag e distance b etw een the E arth and the Sun). In ea ch case, there is an additional 2 E arth-mass planet (whic h we assume is not observed to tra nsit) following a slightly eccentric orbit ( e 2 = 0 . 1) and an initia l mean star-plane t sepa ration of: 0.0 7 8AU (upper left), 0.08 0AU (uppe r r ight) , 0.080 A U (middle left), 0.0 82AU (middle rig ht ), 0.084AU (lo wer left) o r 0.086 A U (low er rig ht). These s eparations a re near the lo cation of the 1 :2 mean motion reso nance ( ≃ 0 . 0794AU). Note that the vertical ax is scale changes from r ow to row. Even these small changes in the orbital s eparation r e sult in qualitative changes in short and long-ter m structure o f the transit timing v ariatio ns . Here we show simulated data from full n-b o dy integrations with no data gaps and no measurement uncer tainties. In practice, even Kepler misses some transits (e.g., due to data downlink with Earth, spac e craft a bno rmalities) and the transit timing measurements hav e nois e o f ∼ 10 − 3 − 10 − 2 days, dep ending pr ima rily on the br ightness of the ta r get sta r and the size of the planet. F o r so me planetar y systems (e.g., Kepler - 9 b& c Holman et a l. 2010) or triple s ta r sy stems (Car ter et al. 2011; Slawson et al. 2011; Steffen et al. 2011), the amplitude o f the tr ansit (eclipse) timing v ariations is m uch la rger than Kepler ’s measurement uncertainties. F or other planetar y systems (e.g., Kepler-1 1 b- f Lis sauer et al. 2011b), the ma gnitude of transit timing v ar iations are compara ble to the measurement uncertainties. W e expec t the Bayesian a pproach and our s urroga te mo del, to b e most useful for such sys tems, once a sufficient num b er and timespan of obs e r v atio ns hav e b een co llected. F or many systems with no detectable transit timing v ariations, sta tistical metho ds such as thos e describ ed here will play a n impo rtant ro le in establishing the significa nce of non-detection and the implied upp er limits for the mass of any per tur bing planets (e.g., Steffen and Agol 200 5). F o rd, Mo or head & V eras 23 Figure 2: This figure is similar to Figur e 1, except the the o uter planet has b een mo ved closer to the transiting planet. In each case, there is a 2 Ear th-mass planet (which we assume is no t obs erved to transit) fo llowing a slightly eccent ric orbit ( e 2 = 0 . 1) and an initial mean star-planet separation of: 0.062 A U (upper left), 0.063AU (upper right), 0.06 4 A U (middle left), 0.065AU (middle right) , 0.066AU (lower left) or 0 .067AU (low er r ight). These s eparations ar e near the loca tion of the 2:3 mean motion resonance ( ≃ 0 . 0655 A U). Note that the vertical axis scale changes from row to row. Again, e ven these small changes in the orbital sepa r ation r esult in qualitative changes in s hort and long-term str uc tur e o f the transit timing v ar iations. 24 Ba y esian Surrogate Mo del for Ti me Series Analysis Figure 3: In this figure, we co nsider a simulated da ta set ba sed on the b est-fit orbital parameters for the exoplanet HD 162020 . W e sho w the marginalized p os terior pr o b- ability distr ibutions for P 1 / 2 = 2 /f 1 (red; very narrow distribution) and P 2 = 1 / f 2 (blue; bro ad distribution) from a sur rogate mo del with N f = 2. Ba s ed on the F ourier expansion of the Do ppler signature of a pla net on a Kepler ian or bit, one would exp ect that the sur rogate mo del would yield marg ina lized p os terior probability distr ibutions for f 2 which o verlaps mar ginalized po sterior probability for f 1 / 2. W e find that this is not necessarily the case, even for this simulated data set with high signal-to- noise, uncorrela ted Gaussian measur e ment errors and accura te e s timates of the mea surement uncertainties. F or this calculation, we hav e used the a ctual times of obse r v ations . F o rd, Mo or head & V eras 25 Figure 4: This figure is a nalogous to Figure 3, except that the observ a tion times a re chosen randomly . In this cas e , the po sterior for P 2 = 1 / f 2 (blue; broad distribution) ov erlaps with the p os terior for P 1 / 2 = 2 f 1 (red; narr ow distribution), as exp ected for a Keplerian orbit. Contrast this with Figure 3 which us es actual observ ation times. W e co nclude tha t for realistic Doppler da ta s ets, aliasing due to limited n umber of the unevenly spaced o bserv a tions can result in the marginal p oster ior distributions for f 1 / 2 and f 2 not ov erlapping . This p os es a c onsiderable challenge to tes ting the one-pla ne t mo del ba sed on the dominant frequencies.