The Elliptical Ornstein-Uhlenbeck Process

The Elliptical Ornstein- U hlenbeck Process Adam M. Sykulski 1 , Sofia C. Olhede 2 and Hanna M. Sykulska-Lawrence 3 Abstract W e introduce the elliptical Ornstein-Uhlenbeck (OU) process, which is a generalisation of the well-kno wn uni v ariate OU process to bi variate time series. This process maps out elliptical stochastic oscillations over time in the complex plane, w hich are observed i n many applications of coupled biv ariate ti me series. The appeal of the model is that elliptical oscillations are generated using one simple first order stochastic differen tial equation (SDE ) , whereas alternati ve models r equire more compli- cated v ectorised or higher order SDE representations. The second useful feature is that parameter estimation can be performed semi-parametrically in the frequency domain using the W hittle Likelihood. W e determine properties of the mode l including the conditions for stationarity , and the geometrical structure of the elliptical oscillations. W e demonstrate the utility of the model by measuring periodic and elliptical pro perties of Earth’ s polar motion. K eywords: Oscillations; Complex-v alued; W idely Linear; Whittle Like lihood; Polar Motion 1 Introd uction Complex-valued representations of biv ariate time series are wide ly used in statistics [17, 40, 42], signal processing [ 30, 34], and numero us application disciplines [2, 13, 43]. A key adv antage of the com plex-valued representatio n is that it can be conveniently used to separate structur es in cou pled biv ariate time series that are circular or n oncir cular when v ie wed in the com plex plane. In signal pr ocessing, this d ic h otomy is sometime s referred to as pr oper or imp r op er , when spec ifically describ in g the geome tr y of the secon d-orde r structure of time ser ies [30]. A type o f noncircu larity o f particular interest is that of elliptica l oscillations in a biv ariate time ser ie s trajectory , which ar e observed acro ss numer ous applications inclu d ing oceanog raphy [18], seismolog y [32], and planetary geophy sics [4]. W e introduce a process that can mo del such elliptical oscillations in continu ous time. Specifically , we propo se the elliptical OU pr ocess g i ven by the following first order SDE dz ( t ) = ( − α 1 + iβ 1 ) z ( t ) dt + ( − α 2 + iβ 2 ) z ∗ ( t ) dt + dW ( t ) , (1) where z ( t ) = x ( t ) + iy ( t ) , i ≡ √ − 1 , and z ∗ ( t ) is the complex conju gate of z ( t ) . The param eters { α 1 , β 1 , α 2 , β 2 } are real-valued, and we shall place con straints on the se p arameters wh ic h en sure z ( t ) is stationary in Section 2. In (1), W ( t ) is a com plex W iener process, whose incre ments follow a co mplex norm al distribution su c h that B = { W ( t + δ ) − W ( t ) } / √ δ ∼ C N (0 , σ 2 , r ) , where σ 2 = E( B B ∗ ) defines th e variance of the complex no r mal, and r = E { B 2 } defines the pseudo- variance and is a complex-valued quantity in general [3 0]. If we set α 2 = β 2 = r = 0 in (1) then we recover the complex OU p rocess, in troduced by Arat ´ o et al. [1 ], which is a cir c ular and prop er com plex-valued p rocess. A proper proce ss form ally mean s that the comp lementary covariance defined by r z ( τ ) = E { z ( t ) z ( t + τ ) } is z e ro for all τ , wher e the autocovariance of a com p lex-valued process is defined by s z ( τ ) = E { z ( t ) z ∗ ( t + τ ) } . Setting α 2 = β 2 = r = 0 in (1), and hence r z ( τ ) = E { z ( t ) z ( t + τ ) } = 0 , has the effect o f ensu ring the comp lex OU of [1] maps out stoch astic circu lar oscillations with frequ ency β 1 and d amping α 1 > 0 . The comp lex OU was pr oposed by [1] specifically to study the Ch a ndler wobble—a small oscillatory deviation in th e Earth’ s axis o f rotatio n , but has also be e n used in nu merous other ph ysical application s in cluding p hysical ocean o graph y [ 33], mag netic fields, a n d rea ction-diffusion systems [2]. The p urpose of this pap er is to study equation (1) in the more gen eral case α 2 6 = β 2 6 = r 6 = 0 . In Sectio n 2 we d eriv e proper ties including co nditions fo r station arity , th e an alytical f o rm of the power spectr al den sity , and the geom etrical relation ship between the SDE para m eters and th e properties of the elliptical oscillations (e. g. the ecc e ntricity and or ientation). I n Section 3 we p rovide computation ally-efficient techniq ues for fitting para m eters of our model to sampled time series, either using a fu lly parametric approach , or a sem i-parametric approach wh e n the mode l is misspecified at certain frequ encies. Finally , in Section 4 we demonstrate the applicability of our model by studying the elliptical oscillations contained within Earth’ s p o lar motion, thus extending the earlier analy ses of Arat ´ o et al. [ 1] and Brillinger [ 5] who restricted findin gs to capturing the pr operties of strictly circular oscillation s. 1 Departmen t of Mathemati cs and Statistics, L ancast er Univ ersity , UK (email: a.sykulski @lancaste r .ac.uk) 2 ´ Ecole polytec hnique f ´ ed ´ erale de Lausanne, Switzerla nd 3 Departmen t of Aeronautic al and Astronautica l Engine ering, Uni versit y of Southampton, UK 1 1.1 Relationship to Literatur e The literature on stochastic mode lling of nonc ircular and improp er complex-valued time series h as primarily f ocused o n linear filters of discrete-time processes, see e.g. [2 3, 28, 32]. T o c r eate noncirc ularity/impr o priety these filters take a wid e ly linear form by ap plying autor egressi ve and moving av erage term s to com p lex-valued p r ocesses and their complex co njugates, taking the ge neral form z t = p X j =1 g j z t − j + p X j =1 h j z ∗ t − j + q X j =0 k j ǫ t − j + q X j =0 l j ǫ ∗ t − j , (2) where ǫ t is i.i.d complex-pro per no ise. In this context o u r generating SDE of (1 ) can be in terpreted as the continuous-time analogu e of the AR(1) version of (2) (with p = 1 , q = 0 ) studied in [32]. This is con sistent with OU p rocesses be in g consider ed the continuous-time analo g ue of AR(1) pr ocesses in general. Howev er, as we shall see, the mapping b etween ( 1) an d (2) is non-tr ivial meaning the pro cesses are worth studyin g separately in their own right—as has been shown to gen erally be the case between CARMA (con tin uous-time ARMA) and discrete-time ARMA models [6 , 8]. In th is paper we focu s on continuou s-tim e pr ocesses, as in many ph ysical app lications it is preferable to mod el the ev olution of a tim e ser ies in continu ous time using SDEs, rather than discrete-time filters. This is b ecause SDE repr esentations allow explicit con nections to be made with u nderlying d ynamical equ a tions (see e.g. [2, 38]), and also provid e a mor e ro bust mo delling framework to deal with h igh frequ ency d ata [7]. In the co ntext of com plex-valued time series, continuous-time mod els hav e been conside r ed in [27] who use Kar hunen - Lo ` eve expansions to ge nerate improp er con tinuous-tim e n onstationary time ser ies. Here we focus o n a station ary mo del and go into depth in terms of und erstanding its statistical prop erties, as well as providing technique s f or param eter estimation, and a demonstration of its app licability to a real-world pr oblem. 2 Pr operties of the elliptical OU process 2.1 Pr ocess Realis ations In Fig. 1 we show two realisations o f the elliptical OU pro cess under two dif ferent sets of parameter values, along with the ir empirical perio dogram estimates to the power spectral density fro m sampled ob servations. The time series are gene r ated usin g th e Euler-Maruyama sch e m e. These realisations explain the u se of the term “elliptical” to describ e the process, as sto c hastic elliptical paths are being traced o ut over time. For complex-valued time series, the periodo gram is in gener al asym metric over positive and negati ve frequen cies, as d irections o f spin are separated in co mplex-valued tim e serie s mod elling. Negative freq u encies correspo n d to clockwise oscillations and positive frequen cies correspo nd to anti-cloc k wise oscillations. In each pa nel we overlay the theo retical p ower spec tr al density whose fu nctional form w ill be derived in Section 2.3, wh ere we also include the effect of aliasing from sampling. As d emonstrated in the figure, th e p rocess acco mplishes ge n erating elliptical oscillations u sin g a simple first order m o del. These elliptical oscillations are seen to have differing e ccentricities, orientations, and ra te s of damping in each example. The oscillations create two peak s o f different magnitu de in the power spectral d ensity , located at th e same correspo n ding n egati ve an d positive f r equency . Equation (1) spe cifies the evolution o r dyn amics of z ( t ) . If we try to decom pose z ( t + δ t ) g i ven z ( t ) , then as dW ( t ) is no t predictable , we have th at z ( t + δ t ) − z ( t ) = Z t + δt t { ( − α 1 + i β 1 ) z ( t ′ ) dt ′ + ( − α 2 + i β 2 ) z ∗ ( t ′ ) dt ′ + dW ( t ′ ) } . (3) W e see d irectly f rom this equation that the incremen t is a widely linear tr ansformatio n of z ( t ) to produce z ( t + δ t ) . Starting from the notio n of comp lex ge o metry [20, Ch. 1], we ca n describe the linear vector space mapp ed ou t by co mplex vectors. In th is space the notion of a line has been r e placed by an ellipse. This ellipse can collap se to a lin e o r a circle under special circumstances, if perturb ed b y dW ( t ′ ) , a s we shall sho rtly show . Thus at every time p oint t ′ , a modification is formed by a d ding an ellipse to the current p o sition z ( t ) to get to z ( t + δ t ) . And as z ( t ) is fixed, if we vie w th e pro cess cond itio nally on its starting point, th e n (3) map s out a sequenc e of superim posed ellipses. T o un derstand the geom etry o f this ellipse we consider a deterministic version of (1) an d (3) whe r e dW ( t ) = 0 . Ex pressing this in term s of x ( t ) and y ( t ) we have that dx ( t ) + i dy ( t ) = {− ( α 1 + α 2 ) x ( t ) + ( β 2 − β 1 ) y ( t ) } dt + i { ( β 1 + β 2 ) x ( t ) + ( α 2 − α 1 ) y ( t ) } dt, (4) such th at the parameter α 1 sets the damping of the process in both x ( t ) and y ( t ) if it is greater than zero—as is the case with the regular r eal-valued OU pro cess. The p arameters { α 2 , β 1 , β 2 } set the geometry of the ellipse of the determ inistic m otion as they cause asymmetric intera c tions between x ( t ) and y ( t ) . As discu ssed already , the e llip se becomes a circle if α 2 = β 2 = 0 , an d this can b e clearly seen fr om (4) when the d amping α 1 is set to ze ro. The oth e r extreme is when the ellipse becomes a line which occurs when β 2 1 = α 2 2 + β 2 2 . T o show this is the case, consider (4) wh ere linear motion occurs when x ( t ) = C y ( t ) (where C is 2 -3 -2 -1 0 1 2 3 -30 -20 -10 0 10 20 30 40 50 -3 -2 -1 0 1 2 3 -30 -20 -10 0 10 20 30 40 50 Figure 1: The to p row displays tw o realisations of the elliptical OU process o f (1) with { α 1 , β 1 , α 2 , β 2 , σ 2 , r } = { 0 . 02 , 1 , − 0 . 5 , − 0 . 3 , 2 , 0 . 6 + i } (le f t) and { α 1 , β 1 , α 2 , β 2 , σ 2 , r } = { 0 . 002 , 0 . 5 , 0 . 3 , 0 . 3 , 0 . 15 , − 0 . 09 − 0 . 09 i } (right). The time serie s z ( t ) is in black, and the { x ( t ) , y ( t ) } com ponents are in grey . W e simulate from t = 0 to t = 100 0 but plot f rom t = 9 00 to t = 1 000 . In the bottom row we d isplay empir ical period o grams in black ( on a decibel scale) of th e fu ll series sampled at in teger values of t , and we overlay the theore tica l p ower spectr al density f rom (13) in r ed, with the aliased version in green. some constant) . Setting the damp ing α 1 = 0 again , f rom (4) we th en solve f or the simu ltaneous equations − α 2 = C ( β 2 − β 1 ) and β 1 + β 2 = C α 2 which yields β 2 1 = α 2 2 + β 2 2 for the special case of linear motion . These special cases will b e verified in Section 2 .2 where we form ally derive the eccentricity of the elliptical oscillations of the stochastic pro cess of (1). This geo metric structur e can b e r elated to time delay embedd ing plots [1 6]. In an emb edding plot, ℜ{ z ( t ) } would be plotted against ℑ{ z ( t ) } acr oss time, and this will capture the dynamic s of th e SDE as enc a psulated by the ellip se geo metry . Finally , note th at ( 3) is a continuous-time specification . I t d e monstrates th at increm ents in the p rocess z ( t ) associated with arb itrary increments δ t are arr i ved a t by a widely linear ope r ation with some noisy o ffset. Equatio n (3) also shows that z ( t ) will trace ou t a continuous-tim e trajectory in the p lane, as specified by { x ( t ) , y ( t ) } . 2.2 Pr ocess Properties T o further understand the properties of (1), we first d efine a cir cu lar real-valued biv ariate OU process ( see also [38]) given b y  d ˜ x ( t ) d ˜ y ( t )  =  − α − β β − α   ˜ x ( t ) ˜ y ( t )  dt + A √ 2  dW 1 ( t ) dW 2 ( t )  , (5) where α > 0 ensures stationarity , β ∈ R sets th e frequen cy o f the circular oscillation, and dW 1 ( t ) and dW 2 ( t ) are indep endent real-valued W iener pro cess incremen ts su ch that { W 1 ( t + δ ) − W 1 ( t ) } / √ δ ∼ N (0 , 1) and { W 2 ( t + δ ) − W 2 ( t ) } / √ δ ∼ N (0 , 1) . W e refer the reader to [ 37] for a mo re general overvie w o f multivariate OU proc e sses. Setting ˜ z ( t ) = ˜ x ( t ) + i ˜ y ( t ) recovers th e complex OU process of [1]. In other words, ( 1) and ( 5) are equiv a le n t whe n α 1 = α , β 1 = β , σ 2 = A 2 , and α 2 = β 2 = r = 0 . W e now transfor m (5) to cr eate elliptical oscillations by defining a new process  x ( t ) y ( t )  = QP  ˜ x ( t ) ˜ y ( t )  , Q =  cos ψ − sin ψ sin ψ cos ψ  , P =  1 ρ 0 0 ρ  . (6) The parameter ρ is a stretching p arameter, an d ψ is a rotatio n param eter , which respec ti vely set the eccentricity and orientation of the elliptical oscillations. F or uniq ueness we restric t 0 < ρ ≤ 1 and − π / 2 ≤ ψ ≤ π / 2 . Note th a t P mu st be applied first in ( 6) for Q to h av e an effect. W e can interp r et (6) as a p hysical defo rmation o f the circular process of (5). 3 T able 1: This table provides a m apping between the parameter s of the elliptical OU p rocess of (1) and the biv ariate p rocess of ( 7). The f unction atan 2 is th e fou r quadr ant inv erse tangen t and sgn is the signum f unction. Bi variate SDE to Elliptical OU Elliptical OU to Bi variate SDE α 1 = α α = α 1 β 1 = β 2  ρ 2 + 1 ρ 2  β = sgn( β 1 ) p β 2 1 − α 2 2 − β 2 2 α 2 = β 2  ρ 2 − 1 ρ 2  sin 2 ψ ρ =  | β 1 |− √ α 2 2 + β 2 2 | β 1 | + √ α 2 2 + β 2 2  1 / 4 β 2 = β 2  ρ 2 − 1 ρ 2  cos 2 ψ ψ = sgn( − β 1 ) 2 atan2 ( α 2 , sgn( − β 1 ) β 2 ) σ 2 = A 2 2  ρ 2 + 1 ρ 2  A 2 = σ 2 √ β 2 1 − α 2 2 − β 2 2 | β 1 | W e now express [ x ( t ) y ( t )] T as a self-con tained biv ariate SDE by combin ing (5) a n d (6) such that  dx ( t ) dy ( t )  = QP  Ω P − 1 Q T  x ( t ) y ( t )  dt + A √ 2  dW 1 ( t ) dW 2 ( t )  , (7) where we use that Q − 1 = Q T and where we define Ω =  − α − β β − α  . Equation (7) is a com plicated vectorised expression fo r gen erating elliptical o scillatio n s, w h ich we contrast with the simpler expression given in (1) u sing the co mplex represen tation. Howev er, (7) is useful for und erstanding the dyna m ics of, and placing parameter constraints on (1), as we shall now show . Spe c ifically , we set z ( t ) = x ( t ) + iy ( t ) and show that ( 7) can th en b e written in the form of (1). T o do this we define the relation ship  x ( t ) y ( t )  = 1 2 T  z ( t ) z ∗ ( t )  , T =  1 1 − i i  . (8) By combining (7) an d (8) we th e n ob tain  dz ( t ) dz ∗ ( t )  = 1 2 T H LT  z ( t ) z ∗ ( t )  dt + T H QP A √ 2  dW 1 ( t ) dW 2 ( t )  , (9) where L = QP Ω P − 1 Q T . The elliptical OU SDE is th en obtained fro m exp a nding (9) and taking the to p row , such tha t we obtain dz ( t ) =  − α + i β 2  1 ρ 2 + ρ 2  z ( t ) dt + β 2  1 ρ 2 − ρ 2  (sin 2 ψ − i cos 2 ψ ) z ∗ ( t ) dt + dW ( t ) , (10) where the incremen t process dW ( t ) is defined by σ 2 = A 2 2  1 ρ 2 + ρ 2  , r = A 2 2  1 ρ 2 − ρ 2  e i 2 ψ . By e q uating the paramete r s in (1) and (1 0) we can o btain an exact one-to -one mapping be twe en the p arameter set { α 1 , β 1 , α 2 , β 2 , σ 2 } of the complex SDE of (1), an d the pa rameter set { α, β , ρ, ψ , A 2 } of th e biv ariate SDE of ( 7). T h e mapping in each dire c tio n is giv en in T able 1. The p a r ameter r , which sets the pseudo -variance of the com p lex-valued increme n t process dW ( t ) , is redundan t and should be set as r = − σ 2 β 1 ( β 2 + i α 2 ) , such that the elliptical OU pro c ess is reduced to five free p arameters f r om ma p ping to an ellip tically transfo rmed bivariate OU process. Setting ρ = 1 in ( 10) (an d T able 1) re c overs the three-param eter comp lex OU of [1] and ( 5). In th e m ore general setting, we observe from T able 1 the simp le relationship th at α 1 = α , meaning α 1 sets the d amping rate of the oscillations in (1), and we thus require α 1 > 0 for the elliptica l OU pr ocess to be stationary . T h e parameters { β 1 , α 2 , β 1 } 4 jointly determine { β , ρ, ψ } (the oscillation frequ ency , eccentricity and orientation) a n d we req uire | β 1 | > p α 2 2 + β 2 2 to create a valid mapping betwe e n the two processes. Th e eccen tr icity of the oscillations is given by ε = p 1 − ρ 4 = s 2 p α 2 2 + β 2 2 | β 1 | + p α 2 2 + β 2 2 , such that la rger values of α 2 and β 2 create more e c c entric o scillations. This fo rmally establishes the geometr ic pro perties of the elliptical oscillation s and verifies the results from Section 2 . 1 that the oscillatio n s are circular when α 2 = β 2 = 0 , and co llapse to a line as α 2 2 + β 2 2 approa c h es β 2 1 . In th e next section we derive the p ower spectral de n sity of the elliptical OU p rocess which will p rovide yet fu r ther intuition on the effect o f the different par a meters. Overall, we see that complex-valued modelling provides a much m o re straig h tforward SDE represen tation of elliptica l o scil- lations than biv ariate modelling , a s shown in (1) a nd Fig. 1. Howe ver , mapping to an under p inning biv a riate process, as in (7), allows us to fur ther un derstand th e g e ometry a nd dynamics of the elliptical oscillations, as well as p lace n ecessary parameter constraints. A similar mapping analysis was performed with discrete-time models in [32] b y equating a widely linear autoregressiv e AR(1) pro cess to a correspo nding biv ariate AR(1) proce ss. The m appings between the p arameters are significantly dif feren t here as compar ed with those found in [32] for discrete time. Th ere are two reasons why these mappin gs are so different. First, althou g h an AR(1) proce ss can generally b e in terpreted as a discrete-time analogu e o f a n OU pr ocess, there is no simple transform ation between their sets of param e te r s in the wid ely linear case, as we show in Appen dix A. This is consistent w ith [ 6, 8] who discuss the nontrivial r elationship between samp led CARMA (co n tinuous-tim e ARMA) mo dels an d regular discrete-time ARMA models. Secondly , the elliptical OU of ( 1) has coefficients given in Cartesian for m, whereas the coefficients of the widely linear AR(1) are giv en in polar form (see (1 9) in Appen dix A). These p arameterizatio n s in each case make sense as the mappings betwee n the O U a n d A R pr ocesses are then straightf orward in th e regu la r (non widely linear) case, see (2 0) in Ap pendix A. Howev er, these choices of p arameterization s cause further departur es in the parameter mappings in the wide ly linear case. As a r esult, the condition s for stationar ity , and the g eometrical proper ties of elliptical oscillation s, are en tirely different in th e co ntinuou s-time elliptical OU prop o sed in this p aper, and th e d iscrete-time AR(1) propo sed in [32]. 2.3 The Po wer Spectral Density For station ary co mplex-valued p rocesses the p ower spectral density can in gen eral be d efined f rom th e autocovariance seq u ence of th e pro cess, such that S z ( ω ) = Z s z ( τ ) e − iωτ dτ , s z ( τ ) = E { z ( t ) z ∗ ( t + τ ) } , (11) where the freq uency ω will always be given in ra dians in this paper . T he p ower spe c tral de n sity of the complex OU of [ 1 ] is given by [ 33] S ˜ z ( ω ) = σ 2 α 2 1 + ( ω − β 2 ) 2 = A 2 α 2 + ( ω − β ) 2 , (12) which we have p rovided b oth in terms of the parameterization of (1) ( with α 2 = β 2 = r = 0 ), and of the circu lar biv ariate process of (5). Note that despite being a proper process, the spectral den sity will contain energy at both n egati ve an d positiv e frequen cies, de caying at rate ω − 2 from the peak frequency . The power spectr al density o f the elliptical OU process is given by S z ( ω ) = A 2       1 ρ + ρ  2 α 2 + ( ω − β ) 2 +  1 ρ − ρ  2 α 2 + ( ω + β ) 2      , (13) which is given in terms of the parame te r ization of the elliptical biv ariate process of (6 ). Then to fin d the power spectral density in terms o f (1) one simp ly substitutes using th e transformations in the r ig ht column of T able 1. T he d eriv ation of (13) is p rovided in App endix B. Intuition is gained by examinin g (13). While (1 2) has just one peak in the spectral den sity loc a te d a t ω = β , (13) has two peaks lo cated at ω = ± β . The rate of damping of both peaks is determ ined by α , and th e ratio of magnitu des of the two peak s is determined by ρ . Note that the orientatio n par ameter ψ does not feature in th e p ower spectral density . When (13) is represented using th e p arameters o f ( 1) then we see that α 1 defines the dampin g of the two p e a ks, a n d { β 1 , α 2 , β 2 } together determin e the peak location s and their relative magnitu des. W e have overlaid the power spectra l density of (13) over the perio dogram o f the simulated series in Fig. 1, where we have also included th e aliased spectr a l density giv en by ˜ S z ( ω ) = P ∞ k = −∞ S z ( ω + 2 π k / ∆) which accounts fo r the departu res at high frequency between the p eriodog ram and spe ctral density , ca used by samplin g the time series a t regular intervals of ∆ . 5 T o fully specify the pro perties of the elliptical OU pro cess, we n eed to derive the complementary spectr u m defined by R z ( ω ) = Z r z ( τ ) e − iωτ dτ , r z ( τ ) = E { z ( t ) z ( t + τ ) } . The comp lex OU of [ 1] is a p r oper process and theref o re R z ( ω ) = r z ( τ ) = 0 . The elliptical OU p r ocess has a co mplementar y spectrum g iven by R z ( ω ) = A 2 4  1 ρ 2 − ρ 2   1 α 2 + ( ω − β ) 2 + 1 α 2 + ( ω + β ) 2  e i 2 ψ , (14) which is dep e ndent on ψ , as well as all the o ther parameter s. W e see that as lon g as ρ < 1 th en R z ( ω ) is no n-zero such that r z ( τ ) is also non -zero and the elliptical OU is an improper pro c e ss as intend e d . The derivation of (14) can also b e found in Ap pendix B. W e note that full specification of the power sp e ctral d ensity and complemen tary spe c tr um allows for an exact method of simulating the proc ess at a fixed sampling rate, based on circu lant embed ding and Fourier transfor ms, as an alterna tive to Eu ler-Maruyama, see [35] f or d etails. 3 Parameter estimation The elliptical OU o f (1) is an impro per proce ss, as we have shown. Theref o re to e stima te p arameters u sing a maxim um like- lihood ap p roach fr om an o bserved complex-valued time series, we would need to in vert large matr ic e s containin g both auto co- variance an d complem entary covariance terms, which will be computatio n ally inten si ve f or large samp le sizes. In this section we detail how inference fo r complex-valued time series can be done in the frequ ency dom ain using the Whittle L ikeliho od and computatio nally-efficient Fourier transform s, see [34] for a compr ehensive review . The Whittle likeliho od is a pseu do-max im um likeliho o d appro ach which has been shown, fo r large c la sses o f p rocesses, to conv erge at th e op tim al O (1 / √ n ) rate to the true param eter values as the sam ple size n increa ses [12]. T he classical assumption s made on the process require bou ndedness (from above and below) a n d twice differentiability of the spectr al density with respect to ω (th e frequen cy) an d θ ∈ Θ (the p arameter vector), as well as the true p arameters lying in th e interior of the parameter space. The elliptica l OU pr o cess indeed has a spec tr al density tha t is bou nded from ab ove and below , and is twice differentiable (with ω and θ ∈ Θ ) , therefor e as long as the tru e parameters do not lie on the bo undary (e.g. linear motion when ρ = 0 ) then we can expect Whittle likelihood to p e rform well for the elliptical OU process with large enoug h sam ples. Consider a length - n observed co mplex-valued time series Z = [ Z 1 , . . . , Z n ] where the time series is regularly sampled at intervals den oted by ∆ . T o perf o rm in ference we need two objects, the first is the Discrete Fourier T ransfor m of th e data Z and its co njugate Z ∗ , den oted J C ( ω ) and given by J C ( ω ) =  J Z ( ω ) J Z ∗ ( ω )  = r ∆ N n X t =1  Z t Z ∗ t  e − iωt ∆ , and the second is the spe c tr al matrix of th e mod el family , deno ted S C ( ω ; θ ) , and given by S C ( ω ; θ ) =  S z ( ω ; θ ) R z ( ω ; θ ) R ∗ z ( ω ; θ ) S z ( − ω ; θ )  , where for the e lliptical OU process, S z ( ω ; θ ) an d R z ( ω ; θ ) ar e as defined in (13) and ( 14) respectively . Th e effect of aliasing can be inco rporated by using ˜ S z ( ω ; θ ) = P K k = − K S z ( ω + 2 π k / ∆; θ ) an d ˜ R z ( ω ; θ ) = P K k = − K R z ( ω + 2 π k / ∆; θ ) in p la c e of S z ( ω ; θ ) and R z ( ω ; θ ) respectively , wh ere the co mputation becomes exact a s K → ∞ , but due to the relatively fast ω − 2 decay in fr e quency , setting a value o f K = 10 w as found to have pr actically conver ged. T o obtain param eter estimates we maximise the following pseud o-likelihood over the parameter space θ ∈ Θ , ℓ W ( θ ) = − 1 2 X ω ∈ Ω  log | S C ( ω ; θ ) | + J H C ( ω ) S − 1 C ( ω ; θ ) J C ( ω )  , (15) where H denotes the Hermitian transpose, and Ω is the set of Fourier frequencies given by Ω = 2 π n ∆ ( −⌈ n/ 2 ⌉ + 1 , . . . , − 1 , 0 , 1 , . . . , ⌊ n/ 2 ⌋ ) . (16) The maximisation of ℓ W ( θ ) is pe rformed over the parameter vector θ = { α, β , ρ, ψ , A 2 } , an d the par ameter estimates corre- sponding to (1) are then fo und u sing the right column o f T ab le 1. Th is app roach can be adapted to irregularly spa c ed observations using the techniqu es described in [22]. 6 A semi-pa rametric alternative is to fit a simple r pseudo- likelihood to only the power spectral d e nsity as gi ven by ℓ S ( θ ) = − X ω ∈ Ω  log S z ( ω ; θ ) + I Z ( ω ) S z ( ω ; θ )  , (17) where I Z ( ω ) = | J Z ( ω ) | 2 is the per iodogra m , and the comp lementary spectru m R z ( ω ; θ ) is n ot used in the fit. Th e aliased spectral d ensity ˜ S z ( ω ; θ ) can be u sed in place of S z ( ω ; θ ) . This semi-parame tric ap proach has so me advantages in terms of robustness to mo del misspecification and smaller sample sizes, as we shall sho rtly d iscuss. This parametric fit, howev er, can only be perfo rmed over the p arameter vector θ = { α, β , ρ, A 2 } as ψ is n o t p resent in the power spectr al density of (13). T o estimate ψ we observe from ( 14) that arg { R z ( ω ) } = 2 ψ , from which we can derive the following no n-param etric estimate ˆ ψ = 1 2 [arg { J Z ( ω max ) } + arg { J Z ( − ω max ) } ] , (18) the full deriv ation of which is given in Appendix C, where ω max refers to the locatio n of the peak in the spectral density , wh ich can be appro ximated fr o m the period ogram if u nknown. As w ith (15), the p arameter estimates from (17) and (1 8) can be expressed in the form of the parameters of (1) using the righ t colum n of T able 1. Both likelihoo d pro cedures (15) and (1 7) can b e m ade fu rther semi-pa rametric by only inclu ding a subset o f freq uencies from (16) in the r espectiv e summations (see also [29]). This is useful when the Fourier transform is con taminated by high frequen cy no ise (and high freq uencies sho uld b e exclud ed from the fit), or the ch o sen mo del is known to only b e co rrect in a narrow r ange of frequ encies, perh aps b ecause an aggr egation of effects has been observed. In deed we shall employ such proced u res in Section 4 to separate the annual and Chandler wobb le oscillations of Earth’ s p olar motio n . Other modificatio ns to Whittle likelihood in cluding taper ing and differencing the time series, or deb iasing the estimates to account for spe c tral blurr in g (see [31] for a revie w). W e did n ot find such modifications to b e needed for the elliptical OU process. This is b ecause the p rocess has a relatively small dynamic ra nge, o wing to the ω − 2 decay in (13), such th a t th e re is only a small amount o f spectr a l blurring pre sent in the perio d ogram. 3.1 Simulation Study W e now perf o rm a simulation study to comp are the bias and r oot-mean -square error (RMSE) of the parameter e stima te s of the elliptical OU p rocess over 1,00 0 Monte Carlo rep licated time series of the model pr esented in th e left p anel of Fig u re 1. T hese parameter values were cho sen as they are of similar magn itude to those ob ser ved in ou r application to Earth’ s pola r motion in the next section. W e com p are four Whittle likelihood ap p roaches: the fu ll p arametric likelihood of (15), the marginal semi- parametric likelihoo d o f (17) (combined with th e non- parametric estimate of ψ in (18)), and sem i- parametric versions of each where o nly a nar rowband of 4 9 Fourier fre q uencies located aroun d each o f the positive and negative pe ak freq uencies are used ( ω ∈ ± [0 . 7 25 , 0 . 89 7] ). The sample size is set to match the app lication in the next section with n = 175 9 . There f ore the semi-parame tric approa c h uses only 98 / 175 9 ≈ 5% of the Fourier fre quencies av ailable in the Fourier transf o rm, b ut a t the frequen cies where most infor mation about the process is con ta in ed. For the first two metho d s we use the appro ximated aliased spectral d ensity with K = 10 . The results are shown in T able 2. W e d isplay bias and RMSE (relative to the tr ue value, expre ssed as a percentag e) fo r each of the p arameters of the ellip tical OU pr ocess: { α 1 , β 1 , α 2 , β 2 , σ 2 } . T h e two appro a ches, (15) and (17), perfor m very similarly (RMSE in the range of 0.5 - 20%), thoug h th e full likelihood of (15) d oes slightly better in estimating { β 1 , α 2 , β 2 } , wh ich are the parameters that define the shape of the ellipse—this is as expected as more inform ation conte n t in the com plementary spectrum has been u sed to fit the par ameters. Howe ver , the m ost challengin g param eters to estimate are the dam ping and amplitu d e { α 1 , σ 2 } ( RMSE 5- 2 0%) and b oth app roaches perf orm similar ly h ere. As we redu ce the ran ge of frequen cies u sed (third and fourth r ows), th e n both app r oaches have slightly worse fits as expected, though the increase in RMSE is mainly observed in the amp litude param eter { σ 2 } (RMSE 20-25 %) as this is the on ly par a meter that lives at a ll fr equencies—inf ormation abo ut { α 1 , β 1 , α 2 , β 2 } (which define the dam p ing, freque ncy , orien tation and e c centricity of the oscillations) liv e at frequ encies in an d around the peak f requency , which is why th e semi-pa r ametric fits generally perfor m very well. In Figure 2 we display the kerne l density estimate plo ts of the deviation of the pa rameter estimates fo r each parameter (relative to its true value) using (17) and ω ∈ [ − π , π ] ( i.e. the second row o f T able 2). The parame te r estimates ar e app roximately Gaussian, and often within 1 0% of the tr u e par ameter value (excep t for α 1 ) , sugg esting th e infer ence technique is robust, at least for the parame te r s a nd sample size con sidered. The m otiv ation be h ind pro posing two estimation tech niques, (15) and (1 7), is because a) two distinct metho ds provid e a validation to ol for num erical optimisers that may sometime s co nverge differently , thu s flagging parameter estimates which may be stuck at bo undaries or lo cal optim a ; an d b) we found that bo th methods will (inevitably) breakdown and struggle to identify all parame te r s when focu sing on a pa r ticularly n arrowband signal. I n g eneral we have foun d optim ising using (17) to b e mo re robust in suc h settings as on ly four p arameters are fitted rather than five. 7 T able 2: Bias an d RMSE (expr essed as a per c e ntage of the true parameter value) for 4 different Whittle inference method s o f the elliptical OU process, co mputed u sing 1,000 Monte Carlo replicated time series of length n = 17 59 α 1 = 0 . 02 β 1 = 1 α 2 = − 0 . 5 α 2 = − 0 . 3 σ 2 = 2 Method Frequencie s Bias RMSE Bias RMSE Bias RMSE Bias RMSE Bias RMSE eq. (1 5) ω ∈ [ − π , π ] 5.4 0 19.17 0.01 0.52 -0.04 1 . 00 0.01 1.42 3.88 5.55 eq. (1 5) ω ∈ [ − π , π ] 5.1 5 19.09 -0 .08 1.42 0.25 4.15 0.41 4.98 3.86 5.92 eq. (1 7) ω ∈ ± [0 . 725 , 0 . 897] 18 .01 3 0.02 0.01 0.61 -0.04 1.1 3 -0.02 1.5 5 19.02 24 .31 eq. (1 7) ω ∈ ± [0 . 725 , 0 . 897] 4.68 26 .66 0.08 0.70 -0.24 1.82 -0 .09 3.44 3 .69 22.6 1 -50 -40 -30 -20 -10 0 10 20 30 40 50 percentage of true value 0 0.05 0.1 0.15 0.2 0.25 0.3 1 1 2 2 2 Figure 2: Kernel den sity estimates o f the deviation of parameter estima te s fr om T able 2 using (17) and ω ∈ [ − π , π ] . Deviations expressed as a percen ta g e of th e true value. 3.2 Parameter Standard Err ors Parameter standar d erro rs and confiden ce intervals can be ob tained from fitting a single time series using bo otstrap resamplin g proced u res developed in [11]. Th is simple pro cedure, developed fo r Wh ittle likelihood estimates, creates bootstrappe d peri- odogr ams by multiplyin g a spectr a l density estimate b y standard independ ent expon e ntially d istributed r a ndom variables at each Fourier freq uency , and then re-estimating parameters with each bootstrap p ed periodo gram to o btain par ameter stand ard errors and confidence intervals. In T ab le 3 we rep ort the param eter standard error s from the Monte Carlo simula tio n of T able 2 f rom (17) and ω ∈ [ − π , π ] . W e also provide the average b ootstrapped parameter standard er rors using two spe c tral d ensity estimators: the raw p eriodog ram and a smooth ed per iodogr am using the Ep ane ˇ cnikov kernel (bandwidth = 0.07 radians). Th e n umber of bootstrap replicates per time series is set to 100. Th e results are displayed in T able 3. Both bootstrap approache s p rovide g ood estimates of the standard errors. In the app lication section we will use the period ogram- based estimator as it is mor e conservative (caused by the highe r variability in the period ogram yielding larger boo tstrap variability). T able 3: Parameter standard errors (a s a percentage of the true value), u sing the same simulation setup as T ab le 2 T e chnique α 1 β 1 α 2 β 2 σ 2 Monte Carlo 17.64 1.36 3.89 4.82 4.04 Bootstrap: Perio dogram 1 7 .42 2.6 9 7. 77 7 .79 3.84 Bootstrap: Epane ˇ cn ikov 11.78 1 .99 5.68 5.69 2.81 4 Earth’ s polar motion Polar motion measures the d eviation of Ear th’ s r otational axis re lative to its crust. I n the lef t pane l o f Fig. 3 we plo t E arth’ s po lar motion from 18 45 to 2021 in o r thogon al x and y direction s, as m easured in milliarcsecon ds (m as), wh ere 100m as cor responds to 8 a deviation of app r oximately 3 metres a t the Earth’ s su r face. Th is data is publicly av ailable fro m the Internation al Earth Rotatio n and Reference Systems Service Earth Or ien tation p roducts 4 . W e o bserve a slow drift in the time series, especially in the y ax is, coupled with clear oscillatory motion. W e are motiv ated to stu d y this dataset in p articular because [1] a lso studied Earth’ s polar motion wh en prop osing the co m plex OU process. Here we can make use o f over 5 0 ye ars’ worth o f ne w d ata to provide u pdated parameter estimates, and test for th e p resence of ellipticity using our elliptical OU p rocess of (1). -5 -4 -3 -2 -1 0 1 2 3 4 5 cycles/year 0 10 20 30 40 50 60 70 80 dB Figure 3: (Left) Earth ’ s polar mo tio n from Dec 1845 to Oct 2 021, measured in r egular interv als of 0.1 year s. (Right) The periodo gram of Ear th’ s polar motion of Fig. 3 w h en represente d as a complex-valued tim e series. The red b and o f freq uencies correspo n ds to the ann ual o scillatio n , and the blue band to the Chandle r wobble. In the right panel of Fig. 3 we plot the pe r iodogr a m of the complex-valued time series z ( t ) = x ( t ) + iy ( t ) . W e d etect three clear p eaks in the periodo gram. The largest at frequen cy zero is due to the drif t. The smallest, at (negative) one cycle per year, is the annual oscillation. The third, at approx im ately -0.8 4 cycles per year is the Ch a ndler wobble, discovered by astrono m er Seth Carlo Chan dler in 1891. W e will stud y th e pr operties of both oscillations u sing th e elliptical OU of ( 1). T o do this we canno t simply loo k at the precise values and location s of the pea k s in th e spectral density—we also n e ed to consider f requenc ies in the vicinity of the peaks, such that we can estimate the dam ping parameter α 1 of the oscillations. W e have marked in blu e and re d in the right panel of Fig. 3 (respectively) the frequ encies we will use to model the Chandler and Ea r th wobb le oscillations respectively . Specifically , the Chandler W obb le is considered in the ran ge -0.97 to -0.70 cycles per year, and the annual oscillation in th e r ange -1.03 to -0.9 7 cycles per year . W e have also mar ked the correspon ding positiv e frequen cies, wh ich will con tain some elevated power if th is compon ent o f the time series has elliptical structure. For visualisation, we ban dpass filter the p olar mo tion time series with box car filters and display the resulting tim e series in Fig. 4. Th e left p anel uses th e b lu e freq uencies from Fig. 3 an d there f ore corre sp onds to the Chandler wobble motion, and th e right panel uses the red freq uencies and therefor e corresponds to the annual o scillation. The presen c e of damping ca n be seen in both oscillations, especially the Chand ler wobble, which motiv ated the con stru ction of the co mplex OU b y A r at ´ o et al. in [1]. No clear ellipticity is seen in the Chan d ler wobble oscillations, but there appea r s to be some present in the ann ual oscillation s. W e will stud y this in more detail b y fitting complex and ellip tica l OU p rocesses to this da ta . No te that the filtering proced u re perfor med h e re (and choice o f filter) is purely for visualisation , and does not im pact the parame ter e stimation, as the perio d ogram of th e unfilter e d time series is u sed in the Wh ittle likelihood. First, we consider the Chandler W obble over n egati ve freque ncies on ly and fit the complex OU o f [1], which corresponds to the ellip tical OU of (1) when α 2 = β 2 = r = 0 . W e fit the p arameters usin g the semi-param e tric Whittle likelihood of (17), using only negative frequen cies in th e interval ω ∈ [ − 0 . 97 , − 0 . 7 ] . The fit of the power spectral density of the comp lex OU in (12) to the per iodogra m is displayed in th e left p anel of Fig. 5. Althoug h the period ogram is variable, it lies within the 95% pointwise confidenc e inter vals of the mo delled power sp ectral den sity almost everywhere. Confidence intervals ar e estimated using th e asymptotic exponential d istribution of the periodog ram. The estimated parameters (to 3 sign ificant figures) are α 1 = 0 . 03 89 ( in units o f y ears − 1 ), β 1 = − 0 . 842 ( cycles per year) and σ 2 = 183 . Using the boo tstrap pro c e dure describe d in Section 3.2 we obtain 95% confiden c e inte r vals o f [0 . 0167 , 0 . 102] f or α 1 , [ − 0 . 847 , − 0 . 8 35] for β 1 , and [11 9 , 266] for σ 2 . In this section we increase the number of boo tstrap rep licates to 10 , 000 such that the approxim a ted confiden ce intervals converge to th e resolution provided. Arat ´ o et al. [1], in their 1962 analy sis, foun d α 1 = 0 . 0 6 and β 1 = − 0 . 83 9 , but the 95% c o nfidence range fo r α 1 was 4 www .iers.org/IERS/EN/Data Products/EarthOri entationData/eop 9 Figure 4: (Left) T he Chan dler wobble over time, which has be e n ban dpass filtered fro m Fig. 3 with a boxca r filter using the b lu e frequen cies h ighlighted in the per iodogr am. (Righ t) The annu a l o scillation which has b een band p ass filtered from Fig. 3 with a boxcar filter using the red frequencies highlig hted in the period ogram. -0.95 -0.9 -0.85 -0.8 -0.75 -0.7 cycles/year 20 30 40 50 60 70 80 dB -0.9 -0.8 -0.7 cycles/year 20 30 40 50 60 70 dB 0.7 0.8 0.9 cycles/year 20 30 40 50 60 70 dB Figure 5: Semi-parametr ic Whittle fits of modelled spectra (red) to the observed period ogram (black) of Earth’ s pola r motion. In the left panel we fit th e comp lex OU spectrum of (12) in th e f requency inter val of -0. 96 to -0.73 cycles per year which captures the Chan dler wobble. In th e righ t two pan els we fit the e llip tical OU sp ectrum of (13) in th e freq uency intervals of -0.96 to -0.73 and 0. 73 to 0 .96 cycles p e r year . In all p anels 9 5 % pointwise co nfidence intervals of th e power spectral density ar e in r ed shadin g. found to be [0 . 008 , 0 . 1 3] which is ov erall co nsistent with ou r estimate s and confidenc e intervals, but we find a sligh tly lower damping parame ter af ter u tilising over 50 years’ worth of new d ata. However , info r mation abou t the damp ing par a meter lies over very few f requenc ies, and is therefor e a challenging parameter to estimate, as also observed by [1] and in our simulation studies of Sectio n 3.1. In other literature, Brillinger [5] also uses a c o mplex OU process like Ara t ´ o et al. , but makes some seasonal correction s, an d also fin ds α 1 = 0 . 06 cycles p er year with a 95% con fidence r ange o f [0 . 006 , 0 . 1 14] . More broadly , ther e still remains an active research deba te on th e rate o f dampin g of the Chandle r wobble [39], where a variety o f g eophysical models have b e en employed to measure this, but a mo re de ta iled comparison with this liter ature is beyond the scop e of this pa p er . W e next fit the ellip tical OU pr ocess to the Chand ler wobble over negative a nd positive fr equencies ω ∈ ± [0 . 7 , 0 . 97 ] , using (17). The fits are displa y ed in the right p anel of Fig. 5. Clearly the fit to positive freque n cies is poor, with n o o bservable peak in th e pe r iodogr am at the exp e cted oscillation fr equency of th e Chan dler wobble, and se veral values ly ing ou tside of the 9 5% pointwise co nfidence inter vals of the mo delled spectrum . Fitting pa rameters u sin g (15) which incor p orates the comp le m entary spectrum did not imp rove the fits. Overall, there is in sufficient evidence to support the presence of a po siti ve-frequ e ncy peak in the Chand ler wobble correspon ding to an elliptical oscillato ry motion. Th is is co n sistent with the literatur e where the Cha n dler wobble motio n has been d e scribed as “ qu asi -circular” with a very low eccentricity in the rang e [0 . 1 , 0 . 2 3 ] in [1 5]. In ou r case, the periodo gram o f the time series is too variable and co ntaminated by o ther artefacts at positive freq uencies, so ou r model is unable to dete c t this low eccen tricity in the Chand ler wobble oscillation. Finally , we fit the elliptical OU pro cess to the an nual oscillation over negative and positiv e fr equencies ω ∈ ± [0 . 9 7 , 1 . 03] using ( 17), where we fix the peak frequ ency β to be 1 cycle per year, leaving just three parameter s { α, ρ, A 2 } to be estimated. 10 -1.02 -1.01 -1 -0.99 -0.98 cycles/year 20 30 40 50 60 70 dB 0.98 0.99 1 1.01 1.02 cycles/year 20 30 40 50 60 70 dB Figure 6: Semi-parametric Whittle fits of the elliptical OU spec tral de nsity o f (13) (red) to the obser ved p e riodogr am (black) of Earth’ s po lar m otion. Fits performed over fre q uency intervals o f - 1 .035 to -0.96 5 and 0.96 5 to 1.03 5 cycles per year, thus capturing the annu al oscillation. 95% pointwise con fidence intervals of the power spectral d ensity are in red shading. This simplification was r equired to make th e o ptimisation f easible as our effective sample size is only 44 (11 po siti ve an d 11 negativ e Fourier f r equencies each co n taining an amp litu de and p h ase) d ue to th e narr ow modelling interval. Infere n ce u sing (15) was not possible in th is e xamp le, as para m eters co n verged to boun dary values. The fitted spectra using (1 7) are displayed in Fig. 6, and this time th e mod el is a go o d fit, and the period ogram comf ortably lies withing the 95% confidence interval bands at all m odelled frequencies. The estimated par ameters of the elliptical OU and their 95% confid ence intervals ar e given in T able 4. Due to the narr owband n ature of the annu al signal, the use o f this mo del and the estimated parame ter values should b e inter p reted with some ca u tion, as reflected by the wid e confidence intervals. The orientation parameter ψ was estimated n on-par ametrically using (18) to be ψ = 0 . 12 5 , and the ecc entricity of the annua l oscillation was estimated to be ε = p 1 − ρ 4 = 0 . 7 82 (with a 9 5% confidenc e inter val of [0 .639,0 .915]) . This is in br o ad agreem ent but som ewhat different from [15] who discover a “significantly elliptic ann ual motion ” in the rang e [0 . 26 , 0 . 4 9] . For c o mparison , a simple non -parametric estimate using (see [32]) ˆ ε = 2 p | J Z ( ω max ) J Z ( − ω max ) | | J Z ( ω max ) | + | J Z ( − ω max ) | , yields ˆ ε = 0 . 53 0 . The h igher values of eccentricity we estimate as compared with [15 ] are likely due to the ir ap p roach of time-windowing into small intervals, versus our approach of co nsidering the entire tim e ser ie s as o ne stoch astic process. Ag ain, a more detailed analysis is b eyond the scope of this paper, howe ver our example h ere serves as a simple proof- of-con cept of the potential app lications of our n ovel elliptical OU SDE model. T able 4: Elliptical OU p arameter estimates and 95 % confid e nce in terval limits α 1 β 1 α 2 β 2 σ 2 estimate 0.024 5 -1.11 0.122 0 .476 28.4 lower limit 0.0091 -1.4 4 0. 066 0.2 58 1 1.4 upper lim it 0.0 929 -1.03 0.2 5 7 1.00 74.7 5 Discussion and Conclusion Oscillations are key featur es of natural a nd h u man-mad e p henomen a. Often we observe lin ked oscillations that map out the same periodic phen omenon . For deterministic phenom ena, su c h have been studied in [1 9, 25], and for stoch a stic phenomen a in [33, 32]. Continuou s-time time series that are imprope r are, a s we have shown, challen ging to describe but possess interpretable multidimen sio n al dynamics [19]. The aim of this paper h a s been to introdu ce a stru c tu red f o rm of multiv ariate depend ence so that stoc hastic elliptical trajectories are map p ed out, just like single oscillations can be co nceptualised as mapping o ut c ir cles. Complex-valued mode ls, such as the elliptical OU pr o cess, provide rich struc tu ral inform ation as we c a n recover th e geom etric features of the ellipse directly from the observations and estimated p arameters. 11 Multiv ariate stochastic pr ocesses have been the focus of inten si ve research in the la st decade [3, 9, 10, 14, 24]. There is much advantage to modelling u n derlying geometr y in time series [26], but that viewpoint exactly cor r esponds as to how the u nderlyin g structure in the o bservations e volves over time. Oscillations are natu ral a s a mode llin g starting po int when studying stationary pheno m ena. The multivariate gener alisation of an oscillation is an observed trajectory from an ellip se [19]. Th is puts an em phasis on th e classes of mod els starting from oscillation s, broa dening to par tially obser ved trajecto ries on th e ellipse. A number of question s remain unre so lved. Our g e n eralisation of the OU mod e l is just one examp le of a statistical model of temporal structur e. The d ifferential equation linkage has been discussed f urther f or other app lications includ ing random fields by [21]. W e e nvision th at similar extensions could be done to their m odel classes. This would build o n the non-p a r ametric statistical work of [ 41]. Furthermo re, we can seek similar extensio n s to triv ariate an d multivariate time series, building stochastic an alogues to the determ in istic a pproach es taken in [19]. Finally , inspired by the works of [1, 2, 5] , the a p plicability o f the elliptical OU process has b e en d emonstrated by the analysis of Earth’ s po lar m otion and the Chandler wobble. In the fu ture, such an alysis can b e repeated on other plane ts such as Mars [36], especially as r ic h er d atasets become available ma k ing such studies m o re f easible. The challenges of real data examp les will stress test ou r model, and sh ow us what ne w featur es and geometrical structures require inco r porating into the mod e l fram ew ork. All results, figures and table s in this paper can be exactly rep roduced u sing the MA TLAB co de av ailable for free download at https://gith ub.com/AdamSy kulski/EllipticalOU . A ppend ix A: Relationship between the OU and AR(1) Processes Consider the widely linea r complex auto regressi ve proc e ss of [ 32] given by Z ( t ) = λe iζ Z ( t − 1) + γ e iφ Z ∗ ( t − 1 ) + ǫ t , (19) with noise variance σ 2 AR and pseud o -variance r AR . Let us now con trast this with th e elliptical OU process of (1). In the simple (prop er) c a se of γ = α 2 = β 2 = r = r AR = 0 then a sampled co mplex OU (at in tervals ∆ ) is like a comp lex AR(1 ) whe r e λ = e − α 1 ∆ , (20) σ 2 AR = σ 2 (1 − e − 2 α 1 ∆ ) 2 α 1 , ζ = β 1 , thus providing a simple mapp ing between the processes. These r elationships can b e d erived b y considering an Euler-Maruyama expansion of the OU: z ( t + 1 /x ) ≈ (1 − α 1 /x + i β 1 /x ) z ( t ) + p A 2 /xB , where x is large and B is a draw fro m a N (0 , 1) such that repeating x ∆ times we have z ( t + ∆)(1 − α 1 /x + iβ 1 /x ) x ∆ z ( t ) + p A 2 /x x ∆ − 1 X k =0 (1 − α 1 /x ) k B , and then taking x → ∞ we get the relation ships above. In the general case γ 6 = α 2 6 = β 2 6 = r 6 = r AR 6 = 0 then the Euler-Maruyama exp a nsion becomes z  t + 1 x  ≈  1 − α 1 x + iβ 1 x  z ( t ) −  α 2 x − iβ 2 x  z ∗ ( t ) + r 1 x B , where B is a draw from C N (0 , σ 2 , r ) . Th e n repeating x ∆ times and taking x → ∞ we ob serve that λe iζ = lim x →∞ x ∆ / 2 X k =0  1 − α 1 x + iβ 1 x  x ∆ − 2 k  α 2 x − iβ 2 x  2 k  x ∆ 2 k  , γ e iφ = lim x →∞ x ∆ / 2 X j =1  1 − α 1 x + iβ 1 x  x ∆ − 2 j +1  α 2 x − iβ 2 x  2 j − 1  x ∆ 2 j − 1  , which have no clear analytical solutions, such that we can obser ve the nontrivial mapping between the pro cesses in the widely linear c a se. 12 A ppend ix B : Po wer spectral density deri vation T o d erive the power spectral density of the ellip tical OU pro cess, we start from the power spectral den sity o f th e complex OU in (1 2) and con vert to Cartesian forms using the relation ships given in [34] S ˜ x ( ω ) = 1 4 { S ˜ z ( ω ) + S ˜ z ( − ω ) } + 1 2 R{ R ˜ z ( ω ) } , S ˜ y ( ω ) = 1 4 { S ˜ z ( ω ) + S ˜ z ( − ω ) } − 1 2 R{ R ˜ z ( ω ) } , S ˜ x ˜ y ( ω ) = 1 2 I { R ˜ z ( ω ) } + i 4 { S ˜ z ( ω ) − S ˜ z ( − ω ) } , where S ˜ x ˜ y ( ω ) is the cross-spectral density b etween ˜ x ( t ) and ˜ y ( t ) , a n d R{·} and I {·} denote th e real and imaginar y part respec - ti vely . Substitutin g in ( 1 2), and using that R ˜ z ( ω ) = 0 as the co mplex OU is a proper pro cess, we obtain S ˜ x ( ω ) = A 2 4  1 α 2 + ( ω − β ) 2 + 1 α 2 + ( ω + β ) 2  , (21) S ˜ y ( ω ) = A 2 4  1 α 2 + ( ω − β ) 2 + 1 α 2 + ( ω + β ) 2  , (22) S ˜ x ˜ y ( ω ) = iA 2 4  1 α 2 + ( ω − β ) 2 − 1 α 2 + ( ω + β ) 2  . (23) Note that S ˜ x ( ω ) = S ˜ y ( ω ) . Next we find th e po wer spectral densities of the ellip tically transform ed biv ariate O U p rocess of ( 7). First we note b y expand ing (6 ) that x ( t ) = 1 ρ ˜ x ( t ) cos ψ − ρ ˜ y ( t ) sin ψ , y ( t ) = ρ ˜ y ( t ) cos ψ + 1 ρ ˜ x ( t ) sin ψ . This clarifies the g eometric inter pretation of P and Q in (6). It then follows that S x ( ω ) = cos 2 ψ ρ 2 S ˜ x ( ω ) + ρ 2 sin 2 ψ S ˜ y ( ω ) − cos ψ sin ψ S ˜ x ˜ y ( ω ) − cos ψ sin ψ S ∗ ˜ x ˜ y ( ω ) , (24) S y ( ω ) = sin 2 ψ ρ 2 S ˜ y ( ω ) + ρ 2 cos 2 ψ S ˜ x ( ω ) + cos ψ sin ψ S ˜ x ˜ y ( ω ) + cos ψ sin ψ S ∗ ˜ x ˜ y ( ω ) , (25) S xy ( ω ) = cos ψ sin ψ ρ 2 S ˜ x ( ω ) − ρ 2 cos ψ sin ψ S ˜ y ( ω ) + cos 2 ψ S ˜ x ˜ y ( ω ) − sin 2 ψ S ∗ ˜ x ˜ y ( ω ) , (26) where we have u sed that S ˜ y ˜ x ( ω ) = S ∗ ˜ x ˜ y ( ω ) . Substituting (2 1)–(2 3) into (24)–(26) we obtain S x ( ω ) = A 2 4  cos 2 ψ ρ 2 + ρ 2 sin 2 ψ   1 α 2 + ( ω − β ) 2 + 1 α 2 + ( ω + β ) 2  , (27) S y ( ω ) = A 2 4  sin 2 ψ ρ 2 + ρ 2 cos 2 ψ   1 α 2 + ( ω − β ) 2 + 1 α 2 + ( ω + β ) 2  , (28) S xy ( ω ) = A 2 4  cos ψ sin ψ ρ 2 − ρ 2 cos ψ sin ψ   1 α 2 + ( ω − β ) 2 + 1 α 2 + ( ω + β ) 2  (29) + iA 2 4  1 α 2 + ( ω − β ) 2 − 1 α 2 + ( ω + β ) 2  . W e then convert back to com plex using the relationships given in [34] S z ( ω ) = S x ( ω ) + S y ( ω ) + 2 I { S xy ( ω ) } , (30) R z ( ω ) = S x ( ω ) − S y ( ω ) + 2 i R{ S xy ( ω ) } . (31) Substituting (2 7)–(29) into (30)–(31) we obtain S z ( ω ) = A 2 4  1 ρ 2 + ρ 2   1 α 2 + ( ω − β ) 2 + 1 α 2 + ( ω + β ) 2  + A 2 2  1 α 2 + ( ω − β ) 2 − 1 α 2 + ( ω + β ) 2  , R z ( ω ) = A 2 4  1 ρ 2 − ρ 2  (cos 2 ψ − sin 2 ψ )  1 α 2 + ( ω − β ) 2 + 1 α 2 + ( ω + β ) 2  + iA 2 2  1 ρ 2 − ρ 2  cos ψ sin ψ  1 α 2 + ( ω − β ) 2 + 1 α 2 + ( ω + β ) 2  , which simplify to the form s g iv en in (13) a n d (14). 13 A ppend ix C: Non-parametric est imation of the orientation Here we derive the form o f the no n-param etric estimate given in (18). Fro m (14) we have that arg { R z ( ω ) } = 2 ψ , fr om which we can form the direct estimate 2 ˆ ψ = arg { ˆ R Z ( ω ) } = a rg { J Z ( ω ) J ∗ Z ∗ ( ω ) } = arg { J Z ( ω ) } + a rg { J ∗ Z ∗ ( ω ) } = arg { J Z ( ω ) } + a rg { J Z ( − ω ) } , where we ha ve used the cross-periodo g ram estimator given by ˆ R Z ( ω ) = J Z ( ω ) J ∗ Z ∗ ( ω ) to estimate R z ( ω ) . W e evaluate the above expression at ω = ω max where ω max is the known/estimated loca tio n of the peak in the po wer spectr a l density o f the elliptical OU process, which in the a p plication to the ann ual oscillation of Ear th ’ s po la r motion this peak occu rs at ω max = 1 cycle per y ear . Refer ences [1] M . Arat ´ o, A. N. Kolmogorov , and Y . G. Sinai. Evaluation of the param eters of a complex stationary Gauss-Markov process. Doklady Ak a demii Nauk SSSR , 14 6 :747–7 50, 196 2 . [2] S. Baran, C. Sz ´ a k-K ocsis, a nd M. Stehl ´ ık. D-op timal d esigns for complex Ornstein – Uhlenbeck processes. J ourna l of Statistical P la nning and Infer ence , 197 :93–10 6, 2018. [3] M . Barigozzi, H. Cho, an d P . Fry zlewicz. Simultaneous m ultiple chang e-point and f actor analysis fo r high-dimension al time series. Journal of E c onometrics , 206 (1):18 7–225 , 201 8. [4] Y . Barkin and J. Fer randiz. Ellip tica l Chandler po le motio ns of the Earth and Mars. In EGU Gen eral Assembly Confer ence Abstracts , volume 12, p age 2 9 36, 2010 . [5] D . R. Brillinger . An empir ical inv estigation o f the Chandler wobble an d two proposed excitation p rocesses. Bulletin of the Internation al S tatistical Institute , 45(3):41 3 –434, 197 3 . [6] P . J. Brockwell. Continuou s-time ARMA processes. Handbook of statistics , 19:2 4 9–27 6 , 2001. [7] P . J. Bro ckwell, V . Ferrazzano, an d C. Kl ¨ upp elberg. High-fr equency sampling and kernel estimation for co ntinuou s-time moving av erage processes. Journal o f T ime Series Ana lysis , 34(3 ):385– 4 04, 2013. [8] K . Chan and H. T ong . A note on emb edding a discrete p arameter arma mode l in a continuo u s parameter ar ma model. Journal of T ime Series Analysis , 8(3 ):277– 281, 1987. [9] J. Chang, B. Guo, Q. Y ao, et al. Princip al compon ent a n alysis for secon d -order stationary vector time series. The Annals of Statistics , 46 (5):209 4–212 4, 201 8. [10] Z. Che, S. Puru sh otham, K. Cho, D. Son tag, and Y . Liu. Recurrent neural networks f or multivariate time series with m issing values. Scien tific r eports , 8 ( 1):608 5, 201 8. [11] R. Dahlhaus and D. Janas. A frequency domain bootstrap for ratio statistics in time series analysis. The An nals o f Sta tistics , 24(5) :1934– 1 963, 1 996. [12] K. Dzh aparidze and A. Y aglom. Sp ectrum parameter estimation in time ser ies analysis. Developmen ts in statistics , 4:1–9 6, 1983. [13] A. P . Guillau min, A. M. Sykulski, S. C. Olh ede, J. J. Early , and J. M. Lilly . Analysis of non-station ary mod u lated time series with app lications to oceanog raphic surface flow mea su rements. J o urnal of T ime Series Analysis , 38 (5):66 8–710 , 201 7. [14] D. Hallac, S. V are, S. Boyd, an d J. Leskovec. T oeplitz in verse covariance-ba sed clustering of multiv a r iate time series data. In Pr oceedin gs of the 23r d A CM SIGKDD I nternationa l Confer ence on Knowledge Discovery and Da ta Mining , pages 215–2 23. ACM, 2017. [15] J. H ¨ opfner . Chandler and an nual wobbles based on space- geodetic measuremen ts. Journal of Geodyn amics , 36(3):3 69–38 1, 2003. [16] H. Kan tz and T . Schreiber . Nonlinea r time series ana lysis , volume 7. Cambridge un iv ersity press, 2 0 04. 14 [17] M. I. Knight a n d M. A. Nunes. L o ng memory estimatio n for co mplex-valued time serie s. Statistics an d Computing , 29(3) :517–5 3 6, 2 019. [18] J. Lilly and J.-C. Gascard . W av elet rid ge diag nosis o f time-varying ellip tical sign a ls with a p plication to an oceanic eddy . Nonlinear Pr ocesses in Geo p hysics , 13(5 ):467–4 83, 2006. [19] J. M. Lilly and S. C. Olhede. Analysis o f mod ulated m u lti variate o scillations. I EEE T ransaction s on Signa l Pr ocessing , 60(2) :600–6 1 2, 2 011. [20] I. V . Lind ell. Method s for electr omagne tic fi eld ana ly sis . Oxford Un iv ersity Press, O x ford, U K, 1992. [21] F . L indgren , H. Rue, and J. Lindstr ¨ om. An explicit link b etween Gaussian fields an d Gaussian Mar kov ran dom fields: the stochastic partial differential equ ation ap proach. Journal o f the Royal Statistical So c iety: Series B (Sta tistical Method o logy) , 73(4) :423–4 9 8, 2 011. [22] Y . Matsuda and Y . Y ajima. Fourier analy sis of irregularly spaced d ata on R d . Journal of the Roya l Statistical So ciety: Se ries B (S ta tistical Method ology) , 71 (1):191 –217, 200 9. [23] J. Na varro-Moren o . ARMA prediction o f widely linear sy stems b y using the innovations alg o rithm. IEEE T ransactions o n Signal Pr ocessing , 56(7):3 061–3 068, 2 0 08. [24] F . H. Nieto, D. Pena, and D. Saboy ´ a. Commo n seasona lity in multivariate time series. Statistica Sinica , 26(4) :1389– 1410, 2016. [25] S. C. Olhed e. Modulated oscillatio n s in many dimension s. Philosophical T ransaction s of th e Roya l Socie ty A: Mathemati- cal, Ph ysical and Engineering Scien ces , 371(1 984):2 0110551, 2013. [26] S. C. Olhede and A. T . W alden. Local directional denoising. IEEE T ransactions on Sign al Pr ocessing , 53(12) :4725– 4730, 2005. [27] A. Oya , J. Navarro-Moreno, and J. C. Ruiz-Mo lina. W idely linear simu lation of con tin uous-time complex-valued r andom signals. IEEE Sign al Pr o cessing Letters , 18(9):513 –516, 2011 . [28] B. Picinbono an d P . Ch evalier . W idely linear estimatio n with complex data. IEEE T ransactions o n Sig nal Pr ocessing , 43(8) :2030– 2 033, 1 995. [29] P . M. Robinson et al. Gaussian semiparametr ic estimation o f lon g range d ependen ce. The An nals o f statistics , 23(5) :1630– 1661, 19 95. [30] P . J. Sch reier and L. L. Scharf . Statistical signal p r oce ssing of complex-valued d ata: the th eory of impr o per and noncircular signals . Cambridge University Press, 2010. [31] A. M. Syk ulski, S. C. Olh ede, A. P . Guillau min, J. M. Lilly , and J. J. Early . The deb iased Whittle likelihood. Biometrika , 106(2 ):251– 2 66, 20 19. [32] A. M. Syku lski, S. C. Olhed e , and J. M. Lilly . A widely linear complex auto r egressi ve p rocess o f or d er one. IEEE T ransactions on Signal Pr ocessing , 64(23) :6 200–6 210, 201 6 . [33] A. M. Sykulski, S. C. Olhed e, J. M. Lilly , and E. Dan ioux. Lagran gian time series mo dels fo r ocea n surface drifter trajectories. J ournal of the R oyal Statistical Society: Series C (App lied Statistics) , 65(1):29 – 50, 201 6. [34] A. M. Sykulski, S. C. Olh e d e, J. M. Lilly , and J. J. E arly . Freq uency-dom ain stochastic modelin g of stationary bi variate or complex-valued signals. IEEE T ransactio n s on S ignal Pr ocessing , 65(12):3 136–3 151, 20 17. [35] A. M. Sykulski and D. B. Perc i val. Exact simulatio n of non c ir cular o r impr oper complex-valued stationary Gau ssian pro- cesses using circu lant emb edding. In 2016 I EEE 2 6th In ternationa l W o rkshop o n Machine Learning for Sign al Pr ocessing (MLSP) , pag es 1–6 . IEEE, 2 0 16. [36] T . V an Ho olst, V . Dehan t, and P . Defraigne. Cha ndler wobble an d free co re nutatio n for Mar s. Pla n etary and S pace Sc ie n ce , 48(12 -14):11 45–1151, 20 00. [37] P . V atiwutipon g and N. Phewchean. Alternati ve w ay to derive the distribution of the multiv ariate Ornstein–Uhlenb e ck process. Advances in Differ ence Eq uations , 2 0 19(1) :1–7, 201 9. 15 [38] M. V eneziani, A. Griff a, A. M. Reynolds, an d A. J. Marian o. Ocean ic tu rbulence and stoch a stic mo d els fr om subsurface Lagrang ian d ata for the Northwe st Atlantic Ocean. Journal of Ph ysical Ocean ography , 34( 8):1884 –1906, 2004 . [39] J. V ondr ´ a k , C. Ron , and Y . Chapanov . New d etermination o f per iod and qu ality factor of Chandler wobble, con sid e ring geoph y sical excitations. Advances in Space Researc h , 59( 5):139 5–140 7, 2017 . [40] P . W ah lb erg and P . J. Schre ier . Lo cally stationar y harm onizable com plex impro per stochastic pro c e sses. Journal of T ime Series Ana lysis , 32(1 ):33–4 6, 2011. [41] A. W alden. Rotary comp onents, ran dom ellipses and p olarization: A statistical perspective. Philosophical T ransactions of the Ro yal Society A: Mathema tical, Physical a nd En gineering Scien c e s , 371(19 84):20 110554, 2 013. [42] A. W alker and T . S. Rao. Periodogram an alysis fo r com plex-valued time series. De velopmen ts in T ime Series Analysis , pages 1 49–16 3, 19 9 3. [43] T . Xiong, Y . Bao, Z. H u , and R. Chiong. Forecasting interval time series using a f ully complex-valued rb f neural network with dp so and pso algorith ms. Info rmation Scien ces , 305:77 –92, 2015. 16