Oscillating systems with cointegrated phase processes

We present cointegration analysis as a method to infer the network structure of a linearly phase coupled oscillating system. By defining a class of oscillating systems with interacting phases, we derive a data generating process where we can specify the coupling structure of a network that resembles biological processes. In particular we study a network of Winfree oscillators, for which we present a statistical analysis of various simulated networks, where we conclude on the coupling structure: the direction of feedback in the phase processes and proportional coupling strength between individual components of the system. We show that we can correctly classify the network structure for such a system by cointegration analysis, for various types of coupling, including uni-/bi-directional and all-to-all coupling. Finally, we analyze a set of EEG recordings and discuss the current applicability of cointegration analysis in the field of neuroscience.


Introduction
Since the first scientific discovery of two pendulums synchronizing by Christiaan Huygens in the seventeenth century, this naturally occurring phenomenon has now been observed in diverse areas such as fireflies synchronizing their flashing behavior, a theatre audience applauding after a show and also in chemical and biological systems, such as the brain and the heart beats of a mother and her fetus, where coupled oscillators appear, see also Pikovsky et al. (2001). Due to it's pervasive presence, understanding synchronization is of key interest for researchers to understand biological networks, such as the connectivity of the nervous system, circadian rhythms or the cardiovascular system. To a statistician this presents a fascinating challenge of modelling complex behavior in large scale systems and how to infer the data-generating mechanisms. To this day, synchronization is not fully understood, but has been the centre of research for decades as evident in Ermentrout (1985), Kuramoto (1984), Strogatz (1987Strogatz ( , 2000, Taylor and Holmes (1998), Winfree (1967), even the phenomenon of synchronizing pendulums as observed by Huygens, still attracts attention today, see Martens et al. (2013), Oliveira and Melo (2015). Many innovative ideas have been presented since Winfree (1967) began a mathematical treatment of the subject. When Kuramoto (1984) first presented his model of coupled oscillators, this made a huge impact in the field and spawned a new generation of research on synchronization. Kuramotos model is still considered among one of the most significant advancements in the study of synchronization in oscillating systems as acknowledged by Strogatz (2000), and the study of coupled oscillators still attracts a fair interest from researchers Ashwin et al. (2016), Burton et al. (2012), Fernandez and Tsimring (2014), Ly (2014), Ly and Ermentrout (2011).
A long standing problem in neuroscience is to recover the network structure in a coupled system. This could for example be to infer the functional connectivity between units in a network of neurons from multiple extracellularly recorded spike trains, or how traces of EEG signals from different locations on the scalp affect each other, which we will treat in this paper. To the authors knowledge, this challenge is still lacking a sound statistical framework to model and test for interaction in a system, as well as impose statistical hypotheses on the network structure. For this task, cointegration analysis offers a refined statistical toolbox, where detailed information on the connections can be inferred, such as the direction and proportional strength of the coupling. The theory of cointegration was originally conceived by Granger (1981), and has since then also been the subject of intense research, most notably within the field of econometrics. In the monograph by Johansen (1996), the full likelihood theory for linear cointegration models with Gaussian i.i.d. errors is derived, and a framework for estimation and inference on parameters using the quotient test is presented. This well acknowledged framework is popularly termed the Johansen procedure. Even though cointegration analysis has developed from within the field of econometrics, it may potentially be used for different models outside economics, such as biological models in continuous time as we explore here. It has also been applied in climate analysis, see Schmith et al. (2012).
In this paper, we demonstrate how to apply cointegration analysis to a system of linearly phase coupled oscillating processes. To display the applicability of the method, we present a simulation experiment, where we present a statistical analysis of phase coupled systems with varying network structures, including uni-/ bi-directional and all-to-all couplings. We show that we can identify the proportional coupling strengths and directions given by the estimated cointegration matrix parameter. Our work is inspired by Dahlhaus and Neddermeyer (2012), which also introduces cointegration analysis as a statistical toolbox to neuroscientists and new challenges for researchers in cointegration theory. However, in contrast to Dahlhaus and Neddermeyer (2012), we incorporate the fact that we are dealing with continuous systems and also ensure that the cointegration property of the system is well posed as a linear structure. This approach assures that the conclusion on the interaction in the data is accurate in terms of cointegration.
The paper is composed as follows. In Sect. 2 we define a class of phase coupled oscillators, in Sect. 3 we highlight some cointegration theory for the analysis including an extension to discretely observed, continuous time models. In Sect. 4 we present a statistical analysis of linearly phase coupled oscillating systems and in Sect. 5 we analyze EEG recordings from an epileptic subject experiencing a seizure, previously analyzed by Shoeb (2009). We discuss the model and findings, conclude on the research and give an outlook of the future direction of the research in Sect. 6. Technical details are presented in the appendix.
Throughout we use the following notation and conventions: unless explicitly stated otherwise, time t ∈ [0, ∞) is assumed continuous, and the process (x t , y t ) is assumed observed with corresponding polar coordinates (φ t , γ t ) . Here denotes transposition. For a p × r matrix M, with r ≤ p, we denote the orthogonal complement M ⊥ , a p × ( p − r ) matrix such that M ⊥ M = 0 (zero matrix). Also denote by sp(A) the subspace spanned by the columns of a matrix A, and let rank(A) denote the rank of the matrix, i.e., the dimension of sp(A).

Oscillating systems
Studying biological rhythms corresponds to studying systems of periodical processes. Intuitively we define a single oscillator as a continuous time bi-variate process z t = (x t , y t ) ∈ R 2 , t ∈ [0, ∞), such that z t revolve around some arbitrary center. Such a process can be derived from an equivalent process in polar coordinates (φ t , γ t ) , where φ t ∈ R is the phase process and γ t ∈ R is the amplitude process, such that (1) We then define the process z t to be an oscillator if the phase process has a monotonic trend.

Defining a class of coupled oscillators
Definition (1) naturally extends to a system of coupled stochastic oscillators, where we observe p oscillators that interact, i.e., z t ∈ R 2 p . Define a class of oscillators with phase (φ t ∈ R p ) and amplitude (γ t ∈ R p ) processes given by the multivariate stochastic differential equations (SDE) where f, g : R 2 p → R p are real valued vector functions, possibly depending on both φ t , γ t or constant, dW φ t , dW γ t are multivariate standard Wiener processes and Σ φ , Σ γ ∈ R p× p such that Σ i Σ i is a positive semi-definite covariance matrix for i = φ, γ . Assume the properties of (2) and (3) are such that and E[φ kt ] is monotonically increasing as a function of t for each k = 1, . . . , p, (5) where E[·] denotes the mean. Since γ t = (γ 1t , . . . , γ pt ) are interpreted as the amplitudes of the individual oscillators, Eq. (4) is a natural assumption and Eq. (5) ensures that the individual oscillators actually revolve (anti-clockwise) around the center and that they are not "stuck" in some part of the phase space, i.e., their angular velocities are positive. Note that we have defined the phase-trend as positive, corresponding to counter-clockwise rotation in accordance with the standard interpretation of the phase. However, for a negative trending process, one can either look at −φ t or simply interpret rotations as clockwise.
To emphasize the implication of inducing interaction in a system, for the data generating process (DGP) in the x y-plane, we derive a DGP from (2)-(3), see "Appendix 1". Assuming that Hence, with the definitions (2)-(5) we have introduced a general class of coupled oscillators, where the specifications of f and g define the properties of the system, such as interaction. This broad definition of oscillating systems covers among others the Kuramoto model, see example (Sect. 2.5) below and other standard oscillators such as the FitzHugh-Nagumo and the Duffing oscillator. In this paper we will analyze phase coupled oscillators, and therefore we assume that g k (φ t , γ t ) = g k (γ kt ), such that there is no feedback from the phase process φ t into the amplitude and the k'th amplitude is not dependent on the rest. Hence, interaction in the system is solely through f (φ t , γ t ), such that the phase processes are attracted by some interdependent relation.

Linear coupling
The arbitrary function f enables us to choose any transformation of the variables to obtain a coupled system, including unidirectional coupling between phases or periodic forcing of the system if we extend f to depend on t as well, intermittent synchronization dependent on a threshold in process differences, etc. Studying the general case where f (φ t , γ t ) is nonlinear in φ t and γ t is a complex exercise. In this paper we restrict ourselves to models where f is composed of a linear mapping of φ t and a function of γ t , with components, for a real matrix Π ∈ R p× p and constant vector ω = (ω 1 , . . . , ω p ) ∈ R p . With this restriction, the interaction between oscillators is linear in the phase, and the k'th oscillator is only dependent on the intrinsic amplitude γ kt through h(γ kt ). We will refer to such a system as linearly phase coupled.
Although we impose the linear restriction Π on the interaction between phases, we can still model a broad set of coupling structures as we show with examples below. Since the interaction is given by Πφ t , we note that the coupling strength in the system is given as the absolute values of the entries of Π and that row k of Π define how oscillator k depends on the rest. Note also that ω defines the attracting state for the phase relations, see example (Sect. 2.3) below. Normally h(γ kt ) is restricted to a constant, but in Sect. 4 we will relax this and investigate systems where h(γ kt ) is only approximately linear and has a sufficiently low variance. This implies a misspecified model, but as we will show, we can still identify the coupling structure, although inference on h(γ kt ) itself is less meaningful.

Example: Linearly phase coupled system with a degenerate γ t process
Let f be defined as in (7) and assume that γ t is a constant (positive) process such that where ω, μ ∈ R p are constant vectors. For reduced rank matrices Π (2) is a continuous time cointegrated process (see Sect. 3) and f admits a linearly phase coupled system with intrinsic rotating frequencies μ. Note that if Π = 0 then there is no interaction in the system, and the individual oscillators will rotate according to their own μ k > 0, and we refer to the system as independent.
The linear specification Π(φ t − ω) implies that at most one attracting point can exist. As an illustration of this, assume a system composed of two coupled oscillators, with where 0 < α 1 + α 2 < 2. Since ω * = ω 1 − ω 2 define an attracting state of the phase difference φ 1t − φ 2t , then with ω * = 0 the system is attracted towards being in-phase, whereas ω * = π would imply that the system is attracted towards being in anti-phase.
Considering that neither α 1 , α 2 or ω * depend on time, the system cannot switch to a different attracting regime. To illustrate possible coupling structures, consider again the system of two oscillators and assume that ω = 0. Then with α 2 = 0 and α 1 = 0 the coupling between φ 1t , φ 2t is uni-directional φ 2t → φ 1t where the arrow → denote the direction of interaction. Likewise, if α 1 = 0 and α 2 = 0 then φ 1t → φ 2t . However, if both α 1 , α 2 = 0 then φ 2t ↔ φ 1t and the coupling is bi-directional. In general, if φ kt appears in the expression f l (φ t ) for oscillator l = k, then φ kt → φ lt . If the opposite is true, then φ lt → φ kt and if both directions exist, then φ lt ↔ φ kt . For f k (φ t ) = 0 oscillator k is (oneway) independent from the rest, but it can still possibly influence others.
For systems where γ t is a degenerate process, then Σ γ = 0 and g(φ t , γ t ) = 0. With σ φ k = σ k then (6) simplifies to where f k (φ t ) = j Π k j φ jt + μ k . Note that if Π = 0 then (8) is simply a constant trend and hence (9) is a rotating process. One can show that the eigenvalues of the deterministic drift matrix in (9) in this case are complex conjugates, − σ 2 2 ± iμ, where i = √ −1, implying that the solutions to (9) oscillate for μ = 0. The oscillations are damped by the negative real part, but sustained by the noise term.
When γ t is a constant vector process the properties of the system are fully identified by (2). Furthermore, if the noise level of the phases Σ φ is sufficiently small, we can use the Hilbert transform 1 to derive the phase process φ t from observations of either x t or y t . This is a commonly used technique in signal processing and has been applied to oscillating systems as well, see Dahlhaus andNeddermeyer (2012), Pikovsky et al. (2001). For systems where φ t is very noisy, this method is less applicable.

Example: Winfree oscillator
With these definitions (6) becomes This example is taken from Winfree (2001) and extended with noise and phase interaction, and therefore we will refer to (10) as the (noisy) Winfree oscillator. Note that the formulation of dγ kt implies that the amplitude fluctuates around κ k . Due to this, we can for sufficiently small noise Σ γ insist that γ kt ≈ κ k for k = 1, . . . , p and therefore analyze the Winfree oscillator using the cointegration toolbox, assuming a constant γ t in dφ t . In Sect. 4 we analyze the range of noise, Σ γ , where the cointegration analysis still performs well.

Example: Kuramoto model
then (2) is the Kuramoto model extended with a stochastic noise term, for phase coupled oscillators, where K k j denotes the coupling strength between the k'th and j'th oscillators. In the classic version, K k j = K ∀k, j, such that for a certain threshold K c , then with K > K c the oscillators exhibit synchronization. For an arbitrary γ t process we cannot simplify (6), but with a degenerate γ t we obtain the same expression as in (9) with f k (φ t ) as in (11). For the Kuramoto model f is a nonlinear function, hence it is not directly applicable to a standard cointegration analysis where f is assumed linear. To emphasize this fact, consider the special case p = 2, where the Kuramoto model is particularly simple and (11) can be written explicitly as, where β = (1, −1) and (α 1 , α 2 ) = (K 12 , K 21 ). If φ 1t ≈ φ 2t at t = 0 and the values of α 1 , α 2 are large enough, then φ 1t ≈ φ 2t ∀t, such that β φ t ≈ 0 and we can write a crude linear approximation of the sine function: sin(β φ t ) ≈ β φ t , such that This is a coarse, but linear, approximation of the Kuramoto model and we can perform a cointegration analysis assuming this approximation is satisfactory. However, one must be cautious with this approximation. Consider sin(β φ t ), when In this case sin(β φ t ) ≈ π − β φ t , and hence and we see that not only do we add a term with π , but the interaction also reverses sign. Recall that 0 < α 1 + α 2 < 2 which implies a stationary relation in the system in (12), see Sect. 3.2. In (13) this condition is reversed, in the sense that −2 < α 1 + α 2 < 0 will imply stationarity. If 0 < α 1 + α 2 < 2, (13) leads to an explosive system, which is not covered in this paper. Therefore, an essential requirement for an approximation of the Kuramoto model is a regime switching ability of (2). For a model with this property, we propose that cointegration analysis on a piecewise linear approximation of the Kuramoto model does make sense and can lead to correct conclusions regarding the network structure. In this paper we will not deal with non-linear cointegration of oscillating systems, but leave this direction open for future research. For a statistical analysis of nonlinear cointegrated systems of the form α t β , i.e. time varying, or regime switching α coefficients, see Bec and Rahbek (2004) and Kristensen and Rahbek (2013). Note that with a general coupling constant K k j = K , then the simple linear approximation to the Kuramoto model

Cointegration
Cointegration theory was originally developed for discrete time processes, however the ubiquitous use of continuous time models has inspired development of continuous time cointegration theory, see Rahbek (2004, 2001). In order to present cointegration analysis as a framework for phase-processes, we therefore review some background on discrete time processes before entering into continuous time cointegrated models. The first part of this section is based on Johansen (1996) and Ltkepohl (2005).

Integrated process
Assume that φ n is a discrete time vector autoregressive process, where A ∈ R p× p , ε n is a Gaussian white noise and μ ∈ R p is a deterministic term. The characteristic polynomial for (15) is the determinant of I p − Aζ for ζ ∈ C, where I p is the p-dimensional identity matrix. If the roots of the characteristic polynomial are all outside the unit circle, then the initial values of φ n can be given a distribution such that φ n is stationary, see Johansen (1996).
If the characteristic polynomial of (15) contains one or more roots at ζ = 1, then there is no stationary solution of φ n , and we say that the process is integrated. In particular, see Johansen (1996), P = A − I p will have reduced rank r < p and can be written as P = ab with a, b ( p × r ) matrices of rank r . Moreover, the process φ n is integrated of order one, I (1) with r cointegrating relations b φ n under regularity conditions presented in Sect. 3.2. Note that the order of integration is a stochastic property and hence including deterministic terms in a model does not change the order of integration.
In this paper we will only deal with I (1) processes, so when we refer to φ n as integrated, we implicitly mean that φ n is integrated of order 1.

Cointegrated process
Let φ n = φ 1n , . . . , φ pn ∈ R p and rewrite (15) with P = A − I p as As already noted if det(I − Aζ ) = 0 implies |ζ | > 1 then φ n has a stationary representation (as an I (0) process). In particular, P has full rank p and all linear combinations of φ n are stationary. If the ( p × p)-dimensional matrix P has reduced rank r < p then P = ab with a, b, p × r dimensional matrices of rank r . Moreover, the process φ n is integrated of order one, I (1) with r cointegrating stationary relations b φ n provided ρ(I r + b a) < 1 with ρ (·) denoting the spectral radius. This we refer to as the I (1) conditions in the following. Note that if r = 0 the process φ n is I (1) with no cointegration, while if r = p (and ρ(A) < 1) then φ n is I (0), or p stationary linear combinations exist. Under the reduced rank r , the system is written as, with b containing the r cointegration vectors and a the loadings or adjustment coefficients. Note that the entries of a and b are not uniquely identified, since we can use any non-singular transformation to obtain similar results. Rather we identify the subspaces sp(a), sp(b) ∈ R r , that is, the subspaces spanned by the columns of a, b, where we use the normalization of b in order to identify parameters uniquely. Furthermore, let m ⊥ denote the matrix such that sp(m ⊥ ) is orthogonal to sp(m), then a necessary condition for an I (1) process is that |a ⊥ b ⊥ | = 0. For more on estimation and inference in cointegration models, see "Appendix 2". Rahbek (2001, 2004) derive a cointegration theory for continuous time models, and conclude that for a discretely observed process, using conventional methods for discrete time generally apply to inference on continuous time parameters. Consider (2) with f as in (8) and for simplicity ω = 0. This is a p-dimensional Ornstein-Uhlenbeck process. The exact solution is

Continuous time cointegrated models
Note that for the solution (17) to be stationary, then Π must be a full rank matrix, and all eigenvalues must have a strictly negative real part. This implies that if Π is not of full rank, then φ t is necessarily not stationary. Assuming discrete observations of (17) at equidistant timepoints t 1 = 0 < t 2 < · · · < t N = T with timestep δ = t n − t n−1 , the corresponding vector autoregressive process is such that the difference process can be written as where ε ∼ N (0, Ω) and Results (18) and (19) hold in general for multivariate processes. Thus, to obtain an estimate for the continuous time matrix,Π, from the discrete time estimateP, a logarithmic transformation involvingP is required For a univariate process (20) is unique, however this is not the case for a multivariate process, due to the non-uniqueness of the multivariate logarithm. Because of this, we cannot uniquely identifyΠ , even though we have a unique estimateP.
For a continuous time process φ t , however, Rahbek (2001, 2004) conclude that this is cointegrated if and only if the discretely observed process (18) is cointegrated. In this case P is of reduced rank, and can be decomposed P = ab with a, b ∈ R p×r of full rank r ≤ p. However, it also holds that for a non-singular matrix ξ = (β α) −1 exp(δβ α) − I r ∈ R r ×r and matrices α, β ∈ R p×r , such that given weak conditions on the sampling timestep δ (see below), the following relations hold see Rahbek (2001, 2004). Hence, for continuous time cointegrated processes, we can infer on the number of cointegration relations (rank(Π ) = r ) from discrete time observations, and also identify the subspaces spanned by the columns of α and β. Note however that due to the unidentified scaling ξ , we can only identify the subspaces, but not the parameters α, β themselves. They are only unique up to a scaling (ξ ), even though we have imposed the normalization (23) and thus uniquely identified a and b.
In the numerical part, we will refer to estimates of α and β, implicitly referring to the discrete time estimates. In terms of subspaces, there is no difference between the discrete and continuous time, but in order to interpret the continuous time Π matrix, one must translate the discrete estimate to a continuous estimate using (19).
It is important to note that when working with continuous time models, one must be careful with regard to the relation (19) between discrete and continuous time and the sampling timestep δ. Kessler and Rahbek (2004) refer to this issue as the embedding problem, and to ensure that the continuous time model is appropriate, one must check for exp(δΠ ) in (18) that it is non-singular, i.e., | exp(δΠ )| = 0, and that it has no negative eigenvalues. If this is the case and the underlying process is in fact cointegrated, the results above hold.

Likelihood ratio test for rank(Π) = r
Consider discrete observations (φ t 1 , . . . , φ t N ) from the continuous process (17) and denote by H r the hypothesis H r : rank(Π ) ≤ r for r = 0, . . . , p. Then the set of hypotheses H 0 , . . . , H r is nested, and H p correspond to the unrestricted model. The likelihood ratio test (LRT) that compare H r and H p is applied sequentially for r = 0, 1, . . . , p − 1 and continued until H r against H p cannot be rejected, and thus determine the number of cointegrating relations for φ t . The LRT statistic is given by whereλ i are the solutions to the eigenvalue problem (49), see "Appendix 2". The asymptotic distribution of (22) is non-standard and therefore it must be simulated.
Here, to also improve on small-sample performance we perform bootstrap simulations as presented by Cavaliere et al. (2012), in order to determine critical values. Specifically, given the data . . , M are simulated and for each sequence the LRT statistic LRT * (m) is re-computed. The empirical quantiles of {LRT * (m) } M m=1 are then used for testing. With r determined,β is given by the r eigenvectors corresponding toλ i , i = 1, . . . , r and the parameter estimateŝ α,μ,Σ follow by ordinary least squares estimation as outlined in "Appendix 2".

Inference for α and β
Since we identify subspaces for α and β, then a normalization is necessary to identify the parameters uniquely. Ifβ is known, thenα follows by OLS. Hence, if we impose a normalization onβ, we can identify all parameters. A common normalization, see Johansen (1996) where c = (I r , 0 p−r ×r ) is a p × r matrix andβ is any version of the r eigenvectors corresponding to the r largest eigenvalues. This ensures that Extending the idea of normalization to restrictions for α, β, we can impose such under the hypothesis H r . Assume that rank(Π ) = r and that the parameters α ∈ R p×r , β ∈ R p×r , μ ∈ R p and Σ ∈ R p× p are all unrestricted within their corresponding subspaces, except for normalization (23). Possible hypotheses for α, β are linear restrictions as given by where A ∈ R p×m , ψ ∈ R m×r , B ∈ R p×s , ξ ∈ R s×r . The known matrices A and B represent the linear hypotheses and ψ and ξ are parameters to be estimated. It is also possible to combine the hypotheses for α and β and we denote this H α,β .
As an example, assume a system of 3 oscillators φ t = (φ 1t , φ 2t , φ 3t ) with r = 1. If we believe that φ 3t is independent of φ 1t and φ 2t , we can specify the hypothesis such that This restriction imply that φ 1t and φ 2t do not contribute to the dynamics of φ 3t , and hence that the latter is independent.
If we want to investigate a possible 1:1 coupling between φ 1t and φ 2t , we can specify and obtain Note however, that under H β the interaction between φ 1t and φ 2t also influence φ 3t if α 3 = 0. Hence, the system admits the relations φ 1t ↔ φ 2t , φ 1t → φ 3t and φ 2t → φ 3t , where the restriction β 3 = 0 implies that the last two relations are uni-directional. If we believe that φ 1t and φ 2t are bi-directionally coupled, φ 1t ↔ φ 2t , but φ 3t is independent and does not contribute to either φ 1t nor φ 2t , we can phrase this hypothesis as a combination of (24) and (25). This leads to the restricted matrix Other hypotheses, such as equal or proportional coupling strength or l : n coupling, can be specified by appropriately designing the matrices A and B. Thus, a broad variety of linear hypotheses on the parameter Π = αβ can be investigated, notably inference on the coupling directions and the effect of system disequilibrium on individual oscillators. Evaluation of the hypotheses H α , H β , and H α,β all lead to similar likelihood ratio tests. To calculate the test statistic, solve again the eigenvalue problem (49) for the unrestricted model, and dependent on the restrictions A and/or B obtain eigenvalues λ * i for the restricted model. The LRT statistic is then given by where H 0 is a generic substitute for any of H α , H β , H α,β . Each of these statistics has an asymptotic χ 2 distribution with varying degrees of freedom (df), where m and s are the column dimensions of the matrices A and B, respectively. This shows that once rank(Π ) is determined, statistical inference for α and β can be carried out, relatively straightforward. As for the rank determination, an alternative to the χ 2 approximation for inference on α and β is to perform bootstrapping for the test (26), see Boswijk et al. (2016).

Numerical simulations 4.1 General setup
We perform a series of experiments with a system of p = 3 linearly coupled Winfree oscillators such that z t ∈ R 6 and f (φ t ) = αβ φ t . Hence, for i = 1, 2, 3, we have a DGP with Since we examine simulations from the Winfree oscillator, our cointegration model will be misspecified since the amplitude is not deterministic and linear, but rather stochastic and fluctuating. However, since the amplitude of the Winfree oscillator has a relatively steady level (of course this also depends on the noise level), due to the squared multiplicative term in the amplitude process, we can approximate it as a constant. Hence we will do so in terms of analyzing the phase process as a cointegrating system. This also implies in terms of parameter estimation for the phase process, the estimate of the constant μ is a pseudo estimate of the κ parameter for the amplitude process, and hence we will compare the estimates to the true value of κ. For each experiment we simulate 1.000.000 iterations of the oscillator (10) using the Euler-Maruyama scheme with timestep Δt = 0.0002 and then subsample for Δt = 0.1, thus obtaining N = 2000 (equidistant and discrete) observations of z t for t ∈ [0, 200). Subsampling every 5000th values diminishes the discretization error of the simulation scheme.
We use the same initial conditions, 1, 0, 0, 1, −1, 0) , and parameters κ = (0.75, 1, 1 for all the experiments so that the only varying parameter is the coupling structure. Note that the κ parameter for φ 2t is set equal to φ 3t to obtain similar simulated outcomes for some experiments to investigate whether we can distinguish between interaction and independence between these two. We set the cointegration parameters for each experiment individually to impose different coupling structures, and will refer to the relevant model by it's Π k , k = 0, 1, 2, 3 matrix, where k defines the model structure (see Fig. 1).
The discrete time model fitted to the data is specified as where the estimateP is used to obtainΠ through (20). The reported estimate forμ is scaled by the timestep Δt. Note that μ is not time-dependent and hence this model will fit a constant parameter μ to a varying quantity γ t and thus it is misspecified as mentioned above. Model (29) is estimated for all 4 systems of three oscillators and we report the parametersΠ andμ for each system. The latter is compared to κ, which is the level parameter of the γ t process.
In addition to a cointegration analysis we apply the mean phase coherence measure, see Mormann et al. (2000), bilaterally to the wrapped phases (i.e., φ it ∈ [0, 2π) for i = 1, 2, 3) as an initial measure of synchronization between the phases in the system. If R ≈ 1 this implies synchronization (R = 1 means that oscillators are perfectly phase locked). On the contrary, R ≈ 0 implies that the distribution of phase differences is approximately uniformly distributed on the unit circle. Note that the mean phase coherence measure is symmetrical, like correlation, and therefore it cannot reveal uni-directional coupling.
In order to determine the significance of the R measures, we bootstrapped critical values for the hypothesis R = 0. Hence, these values are the same for all experiments and presented along with the measured R values. We compare the resulting value of R to the conclusion of the cointegration analysis. We use the same seed for all experiments so that the outcomes are fully comparable in terms of stochastic noise dW . First we run a simulation with uncoupled oscillators as  Fig. 1 Graphical representation of the four systems, represented by the Π i , i = 0, 1, 2, 3 matrix. The arrows define the direction of interaction, hence φ 2t → φ 1t implies that φ 2t is influencing φ 1t (unidirectional coupling), and φ 2t ↔ φ 1t denotes bi-directional coupling, i.e. φ 1t , φ 2t influence eachother  (31), the Π 1 column displays the uni-directional coupled model (32), the Π 2 column displays the bi-directional coupled model (33) and the Π 3 column displays the fully coupled model (34) a benchmark, and then continue with coupled systems as presented in Fig. 1. Figure 2 display the x-coordinates for t ∈ [100, 150] from a simulation of these four systems. The data analysis is carried out using the free software package R (R Core Team 2015). The source code for simulation and bootstrapping procedures are written in C++ to decrease the runtime, utilizing the interface package Rcpp for R and linear algebra package RcppArmadillo for C++. The source code is available in the package cods as supplementary material.

Independent oscillators
This experiment is used as a reference example. We set so rank(Π 0 ) = 0 and there is no interaction in the system. Simulating the model and unwrapping the phases, we obtain the top-left plot of Fig. 3.
Visual inspection of the plot could lead to the conclusion that φ 2t and φ 3t are coupled, however the mean phase coherence measure R for the phases indicates that this is not the case. The distribution of the mean phase coherence measure is unknown, but can be approximated by bootstrapping for H 0 : R = 0, that is for no synchronization present. 1000 bootstrap samples yield the reported 5% critical values in parentheses above ≈0.17, thus the mean phase coherence measure suggest no synchronization present, which is the case.
Performing now a rank test for the rank of Π 0 in the system, we obtain the first part of Table 1.
The test does not reject the hypothesis H r : r = 0, thus suggesting that there is no cointegration present in the system. This in turn implies that the oscillators are independent in terms of synchronization, in accordance with the DGP for Π 0 , and with the mean phase coherence measure.

Uni-directional coupling
In this experiment we analyze a system with a uni-directional coupling. Let such that rank(Π 1 ) = rank(αβ ) = 1, and we have the stationary relation φ 1t − φ 2t . Since α 2 = α 3 = 0, then φ 2t and φ 3t are acting independently, whereas φ 1t is influenced by φ 2t . Hence, the only coupling is φ 2t → φ 1t . The unwrapped phases for the simulation of model Π 1 are seen in the top-right of Fig. 3. The dashed lines indicate the independent phases from the top-left of Fig. 3, and we see that phases φ 2t , φ 3t are equal to their independent versions, whereas we now clearly find that φ 1t is attracted towards φ 2t due to the coupling structure in the system. Examining the mean phase coherence in Eq. (30) for the system (note that R(φ 2t , φ 3t ) is equal to the value in the previous section), we find indications of some synchronization between the phases φ 1t and φ 2t in the system compared to R(φ 1t , φ 2t ) in the independent model. The value is significant on a 5% level as seen by the reported critical values, whereas for R(φ 1t , φ 3t ) and R(φ 2t , φ 3t ) the reported values are not. However, the mean phase coherence measure does not recognize the uni-directional coupling as is the case here. Thus, it cannot distinguish between φ 1t → φ 2t , φ 1t ← φ 2t and φ 1t ↔ φ 2t .
Results from the rank test are in the second part of Table 1. Here we see that r = rank(Π 1 ) = 0 is clearly rejected, whereas r = 1 cannot be rejected with a p value of 0.568. This indicates the presence of a single cointegration relation, in accordance with the construction of the model.
Fitting the model with r = 1, we obtain the unrestricted MLE regression estimates in Table 2. The cointegration relations are close to their true values (approximately within 1 standard error), and both α 2 and α 3 are statistically insignificant. Moreover, the estimates of β suggests a 1:1 coupling between φ 1 and φ 2 .
Therefore, we perform a likelihood test for reducing the unrestricted model, with restrictions for both α and β H α,β : α = Aψ, with A = (1, 0, 0) β = Bξ, with B = (1, −1, 0) , so that A fix α 2 = α 3 = 0 and B restricts to a 1:1 coupling. This yields the test statistic 3.617 which is χ 2 distributed with 4 degrees of freedom and hence implies a p value of 0.460. Thus, we recover the true uni-directional coupling structure of the simulated phases. The fitted model is presented in the right of Table 2. The conclusion is that we have successfully identified the coupling structure of uni-directional coupled phases in a three dimensional system, with two independent phases, and one dependent. Since φ 3t is completely independent of φ 1t and φ 2t and r = 1 we can discard φ 3t when interpreting the cointegration in the system. Then we can interpret the cointegration parameter α as the coupling strength and β as the coupling scheme, here 1:1. If we had analyzed different data, with a β estimate close toβ = (1, −n, 0) , we could then identify a n:1 coupling between φ 1t and φ 2t . This can be seen from the fact that in this case α k (φ 1t −nφ 2t ) would be a stationary relation, and thus φ 2t would rotate ≈n times slower than φ 1t .

A bi-directional coupling with one independent oscillator
We now look at a system with Hence, rank(Π 2 ) = 1 and again we have 1 stationary relation φ 1t − φ 2t , but now with only φ 3t independent, and a bidirectional coupling φ 1t ↔ φ 2t .
Simulating the Π 2 -model we obtain the bottom-left of Fig. 3. We have included the dashed lines again, as references for the independent system. If we contrast the bottomleft of Fig. 3 with the top-right of Fig. 3, we now find that φ 1t and φ 2t are attracting each other, and hence they are both different from their independent versions. Since |α 1 | = |α 2 |, their coupling strength is equal, and the coupled phases lies roughly in the middle between the independent ones. If we look at the mean phase coherence measure for the pairwise interactions, we find relatively strong evidence of a coupling between the phases φ 1t and φ 2t , the value is higher than in the uni-directional case and it is (again) significant given the bootstrapped critical values. However, again we cannot distinguish between types of coupling structures.
Performing a rank test for cointegration in the system with Π 2 , we see in the third part of Table 1 that H r : r = 0 is clearly rejected, and we find that the rank of Π 2 is estimated to 1 with a p value of 0.707. Hence, we recover the correct dimension of the column space of β, and fitting a model with r = 1 yields the parameters in the left of Table 3.  The only insignificant parameter for the model is α 3 , which is in accordance with the construction of the Π 2 model. Specifying the hypothesis (1, −1, 0) , and performing a likelihood ratio test for the reduction yields a test statistic of 3.340, which follows a χ 2 with 3 degrees of freedom, and result in a p value of 0.342. The fitted model is given in the middle of Table 3. If we instead of H α,β specify implying that α 1 = −α 2 , we obtain a test statistic of 3.880, with 4 degrees of freedom, and a p value of 0.423. Thus, we can also restrict the model to one where the coupling strengths are equal in magnitude. The fitted model is presented in the right part of Table 3. Summing up, in a system of bi-directional coupled oscillators plus one independent, we can identify the correct coupling between them, including identifying the proportionally equal coupling strength between the coupled phases. Again we identify r = 1, and hence we can interpret the cointegration parameters as before, hence α is the coupling strength, and β the interaction, again 1:1 coupling.

Fully coupled system
We specify a system with full interaction between all phases.  The simulated phases are shown in the bottom-right of Fig. 3. Comparing to the dashed (independent) versions, we now find that all phases are different from their independent versions. It appears as if φ 2t , φ 3t dominate the system, since φ 1t is attracted closer to their independent versions than otherwise, but it is also a two-against one (κ 2 = κ 3 = κ 1 ) scheme, and we roughly observe that φ 1t is attracted 2/3 towards φ 2t , φ 3t , whereas φ 2t , φ 3t are attracted 1/3 towards φ 1t . So by the construction of the system, this behavior seems natural. We find that the mean phase coherence measure indicates bilateral synchronization for all phases, and all values are significant. The rank test also gives clear evidence of cointegration and we identify r = 2, as seen in the bottom part of Table 1, where both the hypotheses r = 0 and r = 1 are rejected. Fitting a model with r = 2 yields the left half of Table 4.
The estimated κ's are close to their respective values, whereas some of the α parameters deviate (more than 1 standard error) from their true values. If we inspect the estimatedΠ and compare with the true Π 3 it looks better. The row sums are close to zero as they should be, and the signs are correct. The proportional coupling strengths are off though, especially between φ 1t , φ 3t , but it seems that Π 3 is relatively well estimated considering the identification issues. Recall that we can determine the subspaces sp(α) and sp(β) for continuous time cointegration models, see Kessler and Rahbek (2004), but that we have problems regarding the scaling of Π (see Sect. 3.3). Inspired by the fitted values, we restrict both matrices α and β we find that the test statistic is 1.73, χ 2 distributed with 4 degrees of freedom, and thus a p value of 0.785. Hence, we can reduce the model to one with restrictions that generates the true structure of Π . The estimated model parameters are presented in Table 4, and the correspondingΠ iŝ Concluding on the fully coupled system, we find that we can correctly identify the dimension of the cointegration relations. We can also determine the coupling structure as given by the parameters α and β. However, interpretation in this experiment is more informative in terms ofΠ, since with r ≥ 2, the interpretation of cointegration parameters is not as intuitive as in the case of r = 1. We obtain a Π estimate that is reminiscent of the true matrix, with the true directions of the coupling, and strengths somewhat close to the actual values. Thus, we can interpret the system as fully coupled, in a simplistic (linear) Kuramoto type model.

Strength of coupling and identification of interaction
In this section, we compare the mean phase coherence measure to the cointegration analysis with respect to interactions in the system. More specifically, we look at how strong the coupling constants in Π must be in order for the two methods to conclude correctly on interaction in the system. We reuse the parameter settings (34) from the fully coupled experiment, but use a scaled Π matrix Π → εΠ, for ε ∈ [0, 1], where ε controls the coupling strength. The higher ε is, the stronger the coupling, and hence the attraction between phases. Note that ε = 0 corresponds to the model Π 0 and ε = 1 corresponds to Π 3 . The p values are calculated using bootstrapping as presented by Cavaliere et al. (2012) to obtain an estimate of the asymptotic distribution of the trace test statistics. The aim is to investigate the rank test for varying ε compared to identification of interaction in the system using the mean phase coherence measure. Since low values of ε implies weak interaction, the expectation is that both methods will produce doubtful results in a low value regime. From the previous experiment on the fully coupled oscillators, the mean phase coherence measure produced low values on identifying the interaction of the system, hence we expect that the rank test will outperform for low values of ε.
The experimental setup is 100 repetitions for each value of ε, and in each repetition perform 500 bootstrap samples to estimate the p value for the hypotheses H r : r = 0, 1, 2. Figure 4 presents the median p values for the rank test and median mean phase coherence measures against ε. The top row of the figure shows the p values for H r : r = 0, 1, 2 respectively, and the bottom row shows the mean phase coherence (R) measures for pairs of φ 1t , φ 2t and φ 3t . The dotted lines indicate the p = 0.05 value, under which we reject the hypothesis. For the mean phase coherence measure, the 95% significance level for the hypothesis R = 0 has been determined numerically using bootstrapping and is indicated by the dotted lines. If the R-measure falls below this line, independence cannot be rejected. Note that the conclusion r ≤ 3 means that Π is of full rank and therefore invertible, hence β = I 3 . Correct conclusions in bold Seen in the top row of Fig. 4, at least half the simulations reject H r : r = 0 for ε > 0.12, and at least half the simulations reject H r : r = 1 for ε > 0.11. The test does not reject H r : r = 2 for around 88% of the simulated samples for any values of ε. Thus, for ε > 0.11, we can conclude that there is interaction present in the system, and in most of the simulations we also recognize the true rank(Π ) = 2.
If we turn to the bottom row of Fig. 4, where the mean phase coherence measures are shown, we find that half the simulations does not reject the hypothesis R = 0 for ε < 0.34, 0.36 and 0.35, respectively, for R(φ 1t , φ 2t ),R(φ 1t , φ 3t ) and R(φ 2t , φ 3t ), thus clearly indicating an inferior detection of interaction for small values of ε equivalent to weak couplings.
Concluding on this experiment, we find that the rank test detects interaction in the system already at relatively weak coupling strengths. In contrast to this, the coupling must be significantly stronger for a sound conclusion on interaction in the system when using mean phase coherence as a measure of interaction. Furthermore, when detecting interaction in the system, the rank test is also very capable of identifying the true rank of the system, despite a misspecified model. Higher sample sizes will of course improve the inference results.

Consistency of the rank estimation
To investigate the consistency of the cointegration algorithm, we performed an experiment with 1000 repetitions of simulations for Winfree oscillators, the uni-directional coupling, the bi-directional and the fully coupled systems, respectively, and evaluating the rank test, using the same setup as in Sect. 4.1. Table 5 present the percentages of conclusions regarding hypotheses H r : r = 0, r ≤ 1, 2, 3, for each model. Comparing with critical values at a 5% level, obtained by bootstrapping, see Cavaliere et al. (2012), we find that comparing the percentage of simulations where the test correctly identifies the cointegration rank of 1 for uni-and bi-directional coupling are 76.8 and 69.8%, respectively, at a 5% significance level. For the fully coupled system the percentage is 85.5%, and for an independent system the percentage is 96.2%.
These results show that identification of interaction in a system of coupled oscillators is quite precise, and the rank is underestimated in ≤2.5% of the simulations for any model. In the case of independent or full interaction, the method is very good, whereas for systems with directed interaction, or interaction among some oscillators the frequency of overestimating the rank is ≈20-25%. This discrepancy seems intu- Table 6 Percentage of conclusions on interaction indicated by the rank test and the mean phase coherence measures, at a 5% significance level for a sample size of 2000 itively correct, since for the latter systems the true model is a subset of the model of higher rank. As before higher sample sizes will of course improve the inference results.
In Table 6 we compare, in percentages, the conclusions on interaction in the systems, for each model. The values for the rank test presented here, are the summed values from Table 5 for r = 0. We find that both methods are very adept in identifying interaction in these systems. The results, however, should be held against the previous section, where the rank test outperformed the mean phase coherence measure for weak coupling strength. Also noting the fact, that the mean phase coherence measure cannot account for uni-directional coupling, our overall conclusion is that in terms of identifying interaction in the system, the methods seem to perform equally well for stronger coupling, whereas in explaining the system architecture, a cointegration analysis leaves us with more information on how the network is constructed.

Analysis of EEG data
Electroencephalography (EEG) signals are recordings from electrodes distributed on the scalp of subjects. The recorded brainwave patterns are, among others, used for diagnosing sleep disorders, coma or epilepsy. A study on 22 subjects experiencing epileptic seizures from the Children's Hospital Boston is presented by Shoeb (2009) with the aim of detecting seizures based on multiple hours of recordings for each individual. Figure 5 displays an EEG recording of a single subject during a period that include a seizure identified by Shoeb (2009) between 2996 and 3036 s. The seizure is marked by two red dashed lines in Fig. 5. The labels for the signals refer to the individual electrodes on the scalp. We analyze the four signals FP1-F7, FP1-F3, FP2-F4 and FP2-F8, where FP refer to the frontal lobes and F refer to a row of electrodes placed behind these. Even numbered electrodes are on the right side and odd numbered electrodes are on the left side. Smaller (larger) numberings imply that the electrode is placed closer to (further from) the center of the scalp. Hence FP1-F7, FP1-F3 are measurements from the left side, with F3 placed closer to the center than F7, and likewise for right side signals FP2-F4 and FP2-F8. The electrodes for these four signals mirror each other on the left/right side of the scalp. We analyze the seizure period of 40 s and the 40 s leading up to the seizure, i.e. we analyze the two intervals [2956; 2996] and [2996; 3036] respectively, and refer to these as prior to seizure and during seizure. With a sample frequency of 256 measurements each second there are a total of 10,240 measurements for each of the four signals during the 40 s intervals. For more details on the data, see Shoeb (2009). The objective is to compare two fitted cointegration models with interaction as in Eq. (8) for each period: discretely observed for t = 1, . . . , 10,240 in each of the two intervals. The phase processes of the four signals are estimated using the Hilbert transform (see Sect. 2.3). Figure 6 shows the four signals in the two periods and their corresponding estimated unwrapped phase processes. Hence the offsets are in [0, 2π) for the individual phase processes in each period. If we had not split the measurements at 2996 s, the phases in the bottom right of Fig. 6 would be continuations of the phases in the bottom left. A visual inspection of Fig. 6 shows that when transitioning to the seizure period, the phases change to a slower pace (the slopes decrease). Also, prior to the seizure all four phases are closer with no clear distinction between right side and left side phases. During the seizure, the phases split in two groups: right and left side respectively.  This indicates that the model regime changes when transitioning into the seizure period. Table 7 shows the mean phase coherence measures bilaterally for the 4 phase processes and the average of these. Comparing the columns we find no clear indication of a change in the phase regime when transitioning into the seizure period based on this measure, the average change is only 7.5%. However, the measure does indicate interaction in the system among all phases. Table 8 displays the results of a rank test procedure for the system of the four EEG phase processes.
In accordance with the indications from the mean phase coherence measure, the conclusion is a clear presence of cointegration during both periods. Prior to the seizure the rank test of r ≤ 2 is close to the usual 5% significance level, hence the p value The rank is determined to r = 2 in both periods, although the conclusion is far stronger during the seizure. The significance of the statistics are found using 5000 bootstrap samples prior to the seizure due the border limit case of around 5%, during the seizure the p value is determined from 2000 bootstrap samples here is determined using 5000 bootstrap samples, in contrast to the 2000 bootstrap samples used in the other interval, as the conclusion here is quite clear with a p value ≈0.62. In both cases we choose the rank r = 2 for the system. The fitted models are presented in Table 9 with the model fit prior to the seizure on the left side and the fit during the seizure on the right side. If we first note the estimated μ i 's, these are larger during the seizure and significantly so for FP1-F3 and FP2-F4, implying that these phase processes exhibit significantly higher intrinsic linear trends during the seizure. On the other hand, directly interpreting the cointegration parameters is not clear. Recall that these parameters specify subspaces, in this case within R 4 . We therefore look at the estimatedΠ matrices in Table 10 to compare the models for each period.
Here we can determine an all-to-all coupling during both periods and the estimated cointegration matrices show a clear difference for the two intervals. Prior to the seizure the right side signals FP2-F4 and FP2-F8 are much less influenced by the feedback in the system, whereas during the seizure both experience a much larger feedback from the left side signals FP1-F3 and FP1-F7 respectively. Surprisingly, the FP2-F8 signal does not seem to impose a large influence in the system in either interval. It is also interesting to note the changing signs in the two matrices. The two left side signals exhibit a positive feedback on themselves prior to the seizure, whereas during the seizure they impose a negative feedback both on themselves and the right side signals. This could possibly be part of an explanation of the slight kink seen in the phases around 3015-3020 s halfway through the seizure.
Concluding on this analysis we find, not surprisingly, a fully coupled 4 dimensional system with a clear change in the trends prior to and during the seizure. We find that during the seizure the interaction in the system is much stronger, suggesting the more distinctive phases shown in this interval. Including this temporal effect into a single cointegration model covering the full period by utilizing regime switching cointegration models, would be an interesting pursuit for future work.

Discussion
In this paper we have investigated the use of cointegration analysis to determine coupling structures in linearly phase coupled systems. Using these techniques we can with On the left side is the estimated matrix prior to the seizure, on the right side is the estimated matrix during the seizure a good precision identify the coupling structure as a subspace for this type of model. A standard measure to identify synchronization in the literature is the mean phase coherence measure. Contrary to this standard measure, we can detect uni-directional coupling, and we can construct and test hypotheses on the model in form of linear restrictions in the estimated subspace. Furthermore, comparing the mean phase coherence measure with the cointegration analysis in Sect. 4.6, we found that cointegration detects interaction in a system more robustly and for weaker coupling strength than does the mean phase coherence measure. Combined with the fact that cointegration does not just provide a level of synchronization, but rather the structure of the synchronization mechanism, this technique can be used to infer system structures in a much more detailed manner. Of course this higher level of information comes at a cost, since the mean phase coherence measure is easily implemented for any system, whereas the cointegration analysis is more involved and time consuming. Due to the linear nature of the cointegration theory used, we are not able to cover more complex models, such as the Kuramoto model. Thus, an important extension for future work would be to allow for nonlinear coupling functions. However, the linear structure appears naturally when considering a linearization around some phase-locked state, such as for systems showing synchrony or asynchrony. Another interesting pursuit is to extend the model framework to include nonlinear deterministic trends, such that also models like the FitzHugh-Nagumo or the van der Pol oscillator would be covered. The model considered in this paper was constructed from the starting point of the phase process in the spirit of the Kuramoto model, and noise was added on this level. Another approach would be to start from a biological model or a reduction thereof and introduce the noise on the DGP. This would also lead to non-linearities both in drift and diffusion of the phase process. Finally, high dimensional systems are a major challenge in the area of coupled oscillators, hence it would only be natural to investigate cointegration properties of high dimensional systems. A system of more than two synchronizing oscillators that are nonlinearly phase coupled, facilitate chaotic behavior since phases can then bilaterally attract and repel each other. When the number of oscillators increase, one quickly ends up with intuitive shortcomings. The number of parameters rapidly increase with the dimension of the system, possibly leading to a desirable reduction to a sparse interaction structure. This is a key issue with the cointegration framework, which take into account all individual oscillators, as opposed to a mean-field approach that does not run into the same curse of dimensionality. The quality of the estimators will rapidly decrease with increasing dimension of the parameter space or numerical problems may arise. This problem might be alleviated by imposing a sparse interaction structure through a LASSO L 1 penalization.
Cointegration to identify coupling of oscillators has been attempted before in a neuroscience context by Dahlhaus and Neddermeyer (2012). There, the Kuramoto model is approximated for strongly phase coupled oscillators by setting sin(φ j −φ i ) ≈ φ j −φ i , since the phase differences are assumed to be small. We have used the idea from Dahlhaus and Neddermeyer (2012) of analyzing the unwrapped multivariate phase process. Contrary to Dahlhaus and Neddermeyer (2012), however, we have not linearized the sine function to replicate Kuramoto, since this will cause a discrepancy when the phase difference of two oscillators is closer to π than 0 (or 2π ). To mitigate this problem, we have instead taken the approach of designing a DGP with the properties we are interested in, and which allows for any phase differences. Furthermore, this DGP enables us to specify a cointegration model that comply with data from this DGP. Although it may not fully comply with a biological model, it can point to where necessary flexibility is needed in order to develop more realistic cointegration models for biological processes. A first attempt to analyze EEG signals with cointegration analysis with linear coupling structures has been presented. The results are promising, and reveal a finer dependence structure characterizing states of seizure and non-seizure in epileptic patients, which in this example was not possible from the simple Mean Phase Coherence measure. To fully explore the potential of the cointegration analysis for EEG signals, it would be useful to extend the model and analysis tools to allow for non-linearities and simultaneous treatment of many traces, as well as time varying coupling strengths.
Summing up, by applying cointegration as a technique to the field of coupled oscillators in biology, we open up for a whole new area of applications for this statistical theory. On the other hand, using cointegration methods, biologists can gain new insights into network structures, being able to fit models and carry out statistical hypothesis testing. If the cointegration framework presented in this paper can be extended to include the standard models currently used in the field, cointegration would prove a powerful analysis tool for researchers.

+
The quotation marks in the description in (36) imply that one can intuitively interpret the system in this way, but the system is not mathematically split into these parts, as clearly f and g both depends on φ and γ , and σ φ k enters in the "amplitude" part. Generalizing (36) we find where R t is a block diagonal matrix of 2 × 2 rotation matrices and Q t is a diagonal matrix of amplitude dependent adjustments. The noise is composed of a sum of two state dependent multivariate processes, where the functions a, b define the noise as given in (36). The time index t has been added in (37) to emphasize the time dependency of the matrices R t and Q t .
Define R p× p M i j = N −1 N n=1 ϒ it n ϒ jt n , for i, j = 0, 1 Note that M i j = M ji . Then the estimate of μ given a and b iŝ μ(a, b) = M 02 − ab12 Define the residuals With these preliminary steps, we obtain the profiled likelihood function log L(a, b, Ω) = − 1 2 N log |Ω| − 1 2 N n=1 (R 0t n − ab R 1t n ) Ω −1 (R 0t n − ab R 1t n ), equivalent to the regression equation Equation (43) is estimated as a reduced rank regression, by solving for eigenvalues. Define Then for a fixed b, we obtain estimates for a and Ω by OLS with b R 1t n as the independent variable, and the likelihood is then maximized as Using the Schur complement for the matrix we find |S 00 − S 01 b(b S 11 b) −1 b S 10 | = |S 00 |||b (S 11 − S 10 S −1 00 S 01 )b|/|b S 11 b|.
Equation (48) is maximized for all p × r matrices by solving for the p eigenvalues λ i in |λS 11 − S 10 S −1 00 S 01 | = 0, such that λ i S 11 v i = S 10 S −1 00 S 01 v i , and the p × 1 eigenvectors v i , i = 1, . . . , p are normalized, v j S 11 v i = 1, for i = j 0, for i = j.
Then for a given r ,b ( p × r ) is given by the r eigenvectors, v 1 , . . . , v r , corresponding to the r largest eigenvaluesλ 1 > · · · >λ r , and the maximum value of the likelihood function with thisb is Since (51) holds for r = 0, . . . , p, where for r = 0 we set sp(b) = {0} and for r = p we set sp(b) = R p , we have solved for all possible ranks r once and for all, and we can form the likelihood ratio test