Contents
List of Contributors 9
1 Handbook of Time Series Analysis: Introduction and Overview 13
2 Nonlinear Analysis of Time Series Data
(Henry D. I. Abarbanel and Ulrich Parlitz) 17
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2 Unfolding the Data: Embedding Theorem in Practice . . . . . . . . 18
2.2.1 Choosing T : Average Mutual Information . . . . . . . . . . . 20
2.2.2 Choosing D: False Nearest Neighbors . . . . . . . . . . . . . 24
2.2.3 Interspike Intervals . . . . . . . . . . . . . . . . . . . . . . . . 30
2.3 Where are We? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.4 Lyapunov Exponents: Prediction, Classification, and Chaos . . . . 31
2.5 Predicting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.6 Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.6.1 Modeling Interspike Intervals . . . . . . . . . . . . . . . . . . 40
2.6.2 Modeling the Observed Membrane Voltage Time Series . . 41
2.6.3 ODE Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . 45
2.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3 Local and Cluster Weighted Modeling for Time Series Prediction
(David Engster and Ulrich Parlitz) 51
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.1.1 Time Series Prediction . . . . . . . . . . . . . . . . . . . . . . 52
3.1.2 Cross Prediction . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.1.3 Bias, Variance, Overfitting . . . . . . . . . . . . . . . . . . . . 53
3.2 Local Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.2.1 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.2.2 Local Polynomial Models . . . . . . . . . . . . . . . . . . . . 57
3.2.3 Local Averaging Models . . . . . . . . . . . . . . . . . . . . . 58
3.2.4 Locally Linear Models . . . . . . . . . . . . . . . . . . . . . . 58
3.2.5 Parameters of Local Modeling . . . . . . . . . . . . . . . . . 58
3.2.6 Regularization . . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.2.7 Optimization of Local Models . . . . . . . . . . . . . . . . . 64
3.3 Cluster Weighted Modeling . . . . . . . . . . . . . . . . . . . . . . . 65
3.3.1 The EM Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 67
3.4 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
Handbook of Time Series Analysis. Björn Schelter, Matthias Winterhalder, Jens Timmer
Copyright © 2006 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim
ISBN: ???????????????

2 Contents
3.4.1 Noise Reduction . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.4.2 Signal Through Chaotic Channel . . . . . . . . . . . . . . . . 70
3.4.3 Friction Modeling . . . . . . . . . . . . . . . . . . . . . . . . 72
3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4 Deterministic and Probabilistic Forecasting in Reconstructed State Spaces
(Holger Kantz and Eckehard Olbrich) 79
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.2 Determinism and Embedding . . . . . . . . . . . . . . . . . . . . . . 81
4.3 Stochastic Processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.4 Events and Classification Error . . . . . . . . . . . . . . . . . . . . . 93
4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
5 Dealing with Randomness in Biosignals
(P. Celka, R. Vetter, E. Gysels, and T. Hine) 101
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
5.1.1 Determinism: Does It Exist? . . . . . . . . . . . . . . . . . . . 101
5.1.2 Randomness: An Illusion? . . . . . . . . . . . . . . . . . . . . 102
5.1.3 Randomness and Noise . . . . . . . . . . . . . . . . . . . . . 104
5.2 How Do Biological Systems Cope with or Use Randomness? . . . 105
5.2.1 Uncertainty Principle in Biology . . . . . . . . . . . . . . . . 105
5.2.2 Stochastic Resonance in Biology . . . . . . . . . . . . . . . . 106
5.3 How Do Scientists and Engineers Cope with Randomness
and Noise? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
5.4 A Selection of Coping Approaches . . . . . . . . . . . . . . . . . . . 110
5.4.1 Global State-Space Principal Component Analysis . . . . . 110
5.4.2 Local State-Space Principal Component Analysis . . . . . . 120
5.5 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
5.5.1 Cardiovascular Signals: Observer of the Autonomic
Cardiac Modulation . . . . . . . . . . . . . . . . . . . . . . . 124
5.5.2 Electroencephalogram: Spontaneous EEG and Evoked
Potentials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
5.5.3 Speech Enhancement . . . . . . . . . . . . . . . . . . . . . . . 132
5.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
6 Robust Detail-Preserving Signal Extraction
(Ursula Gather, Roland Fried, and Vivian Lanius) 143
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
6.2 Filters Based on Local Constant Fits . . . . . . . . . . . . . . . . . . 146
6.2.1 Standard Median Filters . . . . . . . . . . . . . . . . . . . . . 146
6.2.2 Modified Order Statistic Filters . . . . . . . . . . . . . . . . . 148

Contents 3
6.2.3 Weighted Median Filters . . . . . . . . . . . . . . . . . . . . . 151
6.3 Filters Based on Local Linear Fits . . . . . . . . . . . . . . . . . . . . 153
6.3.1 Filters Based on Robust Regression . . . . . . . . . . . . . . 153
6.3.2 Modified Repeated Median Filters . . . . . . . . . . . . . . . 155
6.3.3 Weighted Repeated Median Filters . . . . . . . . . . . . . . . 156
6.4 Modifications for Better Preservation of Shifts . . . . . . . . . . . . 157
6.4.1 Linear Median Hybrid Filters . . . . . . . . . . . . . . . . . 157
6.4.2 Repeated Median Hybrid Filters . . . . . . . . . . . . . . . . 159
6.4.3 Level Shift Detection . . . . . . . . . . . . . . . . . . . . . . . 161
6.4.4 Impulse Detection . . . . . . . . . . . . . . . . . . . . . . . . 163
6.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165
7 Coupled Oscillators Approach in Analysis of Bivariate Data
(Michael Rosenblum, Laura Cimponeriu, and Arkady Pikovsky) 171
7.1 Bivariate Data Analysis: Model-Based Versus Nonmodel-Based
Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171
7.1.1 Coupled Oscillators: Main Effects . . . . . . . . . . . . . . . 173
7.1.2 Weakly Coupled Oscillators: Phase Dynamics Description . 175
7.1.3 Estimation of Phases from Data . . . . . . . . . . . . . . . . 176
7.1.4 Example: Cardiorespiratory Interaction in a Healthy Baby . 178
7.2 Reconstruction of Phase Dynamics from Data . . . . . . . . . . . . 179
7.3 Characterization of Coupling from Data . . . . . . . . . . . . . . . . 183
7.3.1 Interaction Strength . . . . . . . . . . . . . . . . . . . . . . . 183
7.3.2 Directionality of Coupling . . . . . . . . . . . . . . . . . . . . 185
7.3.3 Delay in Coupling from Data . . . . . . . . . . . . . . . . . . 187
7.4 Conclusions and Discussion . . . . . . . . . . . . . . . . . . . . . . . 189
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190
8 Nonlinear Dynamical Models from Chaotic Time Series:
Methods and Applications
(Dmitry A. Smirnov and B. P. Bezruchko) 193
8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193
8.2 Scheme of the Modeling Process . . . . . . . . . . . . . . . . . . . . 194
8.3 “White Box” Problems . . . . . . . . . . . . . . . . . . . . . . . . . . 196
8.3.1 Parameter Estimates and Their Accuracy . . . . . . . . . . . 196
8.3.2 Hidden Variables . . . . . . . . . . . . . . . . . . . . . . . . . 200
8.3.3 What Do We Get from Successful and Unsuccessful
Modeling Attempts? . . . . . . . . . . . . . . . . . . . . . . . 202
8.4 “Gray Box” Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . 203
8.4.1 Approximation and “Overlearning” Problem . . . . . . . . 203
8.4.2 Model Structure Selection . . . . . . . . . . . . . . . . . . . . 205
8.4.3 Reconstruction of Regularly Driven Systems . . . . . . . . . 206
8.5 “Black Box” Problems . . . . . . . . . . . . . . . . . . . . . . . . . . 207

4 Contents
8.5.1 Universal Structures of Model Equations . . . . . . . . . . . 207
8.5.2 Choice of Dynamical Variables . . . . . . . . . . . . . . . . . 210
8.6 Applications of Empirical Models . . . . . . . . . . . . . . . . . . . 211
8.6.1 Method to Reveal Weak Directional Coupling
Between Oscillatory Systems from Short Time Series . . . . 212
8.6.2 Application to Climatic Data . . . . . . . . . . . . . . . . . . 213
8.6.3 Application to Electroencephalogram Data . . . . . . . . . . 215
8.6.4 Other Applications . . . . . . . . . . . . . . . . . . . . . . . . 217
8.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 218
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 219
9 Data-Driven Analysis of Nonstationary Brain Signals
(Mario Chavez, Claude Adam, Stefano Boccaletti and Jacques Martinerie) 225
9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 225
9.1.1 EMD-Related Work . . . . . . . . . . . . . . . . . . . . . . . . 226
9.2 Intrinsic Time-Scale Decomposition . . . . . . . . . . . . . . . . . . 227
9.2.1 EMD and Instantaneous Phase Estimation . . . . . . . . . . 228
9.2.2 Drawbacks of the EMD . . . . . . . . . . . . . . . . . . . . . 230
9.3 Intrinsic Time Scales of Forced Systems . . . . . . . . . . . . . . . . 231
9.4 Intrinsic Time Scales of Coupled Systems . . . . . . . . . . . . . . . 232
9.5 Intrinsic Time Scales of Epileptic Signals . . . . . . . . . . . . . . . 234
9.5.1 Intracerebral Activities . . . . . . . . . . . . . . . . . . . . . . 234
9.5.2 Magnetoencephalographic Data . . . . . . . . . . . . . . . . 235
9.6 Time-Scale Synchronization of SEEG Data . . . . . . . . . . . . . . 237
9.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 238
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 239
10 Synchronization Analysis and Recurrence in Complex Systems
(M. C. Romano, M. Thiel, J. Kurths, M. Rolfs, R. Engbert, and R. Kliegl) 243
10.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243
10.2 Phase Synchronization by Means of Recurrences . . . . . . . . . . . 245
10.2.1 Examples of Application . . . . . . . . . . . . . . . . . . . . 250
10.2.2 Influence of Noise . . . . . . . . . . . . . . . . . . . . . . . . 253
10.3 Generalized Synchronization and Recurrence . . . . . . . . . . . . . 256
10.3.1 Examples of Application . . . . . . . . . . . . . . . . . . . . 258
10.4 Transitions to Synchronization . . . . . . . . . . . . . . . . . . . . . 261
10.5 Twin Surrogates to Test for PS . . . . . . . . . . . . . . . . . . . . . 264
10.6 Application to Fixational Eye Movements . . . . . . . . . . . . . . . 268
10.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 271
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 272
11 Detecting Coupling in the Presence of Noise and Nonlinearity
(Theoden I. Netoff, Thomas L. Carroll, Louis M. Pecora, and Steven J. Schiff ) 277
11.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 277

Contents 5
11.2 Methods of Detecting Coupling . . . . . . . . . . . . . . . . . . . . . 278
11.2.1 Cross-Correlation . . . . . . . . . . . . . . . . . . . . . . . . . 278
11.2.2 Mutual Information . . . . . . . . . . . . . . . . . . . . . . . 279
11.2.3 Mutual Information in Two Dimensions . . . . . . . . . . . 280
11.2.4 Phase Correlation . . . . . . . . . . . . . . . . . . . . . . . . . 280
11.2.5 Continuity Measure . . . . . . . . . . . . . . . . . . . . . . . 281
11.3 Linear and Nonlinear Systems . . . . . . . . . . . . . . . . . . . . . 282
11.3.1 Gaussian Distributed White Noise . . . . . . . . . . . . . . . 282
11.3.2 Autoregressive Model . . . . . . . . . . . . . . . . . . . . . . 282
11.3.3 Hénon Map . . . . . . . . . . . . . . . . . . . . . . . . . . . . 284
11.3.4 Rössler Attractor . . . . . . . . . . . . . . . . . . . . . . . . . 284
11.3.5 Circuit Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 285
11.4 Uncoupled Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 285
11.4.1 Correlation Between Gaussian Distributed Random Data Sets286
11.4.2 Correlation Between Uncoupled AR Models . . . . . . . . . 286
11.4.3 Correlation Between Uncoupled Hénon Maps . . . . . . . . 287
11.4.4 Correlation Between Uncoupled Rössler Attractors . . . . . 287
11.4.5 Uncoupled Electrical Systems . . . . . . . . . . . . . . . . . . 287
11.5 Weakly Coupled Systems . . . . . . . . . . . . . . . . . . . . . . . . 289
11.5.1 Coupled AR Models . . . . . . . . . . . . . . . . . . . . . . . 289
11.5.2 Coupled Hénon Maps . . . . . . . . . . . . . . . . . . . . . . 289
11.5.3 Weakly Coupled Rössler Attractors . . . . . . . . . . . . . . 289
11.5.4 Experimental Electrical Nonlinear Coupled Circuit . . . . . 290
11.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 290
11.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 292
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293
12 Linear Models for Mutivariate Time Series
(Manfred Deistler) 295
12.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 295
12.2 Stationary Processes and Linear Systems . . . . . . . . . . . . . . . 296
12.3 Multivariable State Space and ARMA(X) Models . . . . . . . . . . 300
12.3.1 State Space and ARMA(X) Systems . . . . . . . . . . . . . . 301
12.3.2 Realization of State Space and ARMA Systems . . . . . . . 303
12.3.3 Parametrization and Semi-Nonparametric Identification . . 305
12.3.4 CCA-Subspace Estimators . . . . . . . . . . . . . . . . . . . . 307
12.3.5 Maximum Likelihood Estimation Using Data Driven Local
Coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . 309
12.4 Factor Models for Time Series . . . . . . . . . . . . . . . . . . . . . . 311
12.4.1 Principal Component Analysis . . . . . . . . . . . . . . . . . 312
12.4.2 Factor Models with Idiosyncratic Noise . . . . . . . . . . . . 313
12.4.3 Generalized Linear Dynamic Factor Models . . . . . . . . . 315
12.5 Summary and Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . 316
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 318

6 Contents
13 Spatio-Temporal Modeling for Biosurveillance
(David S. Stoffer and Myron J. Katzoff ) 321
13.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 321
13.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 322
13.3 The State Space Model . . . . . . . . . . . . . . . . . . . . . . . . . . 324
13.4 Spatially Constrained Models . . . . . . . . . . . . . . . . . . . . . . 328
13.5 Data Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 332
13.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 345
14 Graphical Modeling of Dynamic Relationships
in Multivariate Time Series
(Michael Eichler) 347
14.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 347
14.2 Granger Causality in Multivariate Time Series . . . . . . . . . . . . 349
14.2.1 Granger Causality and Vector Autoregressions . . . . . . . . 349
14.2.2 Granger Causality in the Frequency Domain . . . . . . . . . 352
14.2.3 Bivariate Granger Causality . . . . . . . . . . . . . . . . . . . 353
14.3 Graphical Representations of Granger Causality . . . . . . . . . . . 354
14.3.1 Path Diagrams for Multivariate Time Series . . . . . . . . . 354
14.3.2 Bivariate Granger Causality Graphs . . . . . . . . . . . . . . 356
14.4 Markov Interpretation of Path Diagrams . . . . . . . . . . . . . . . 358
14.4.1 Separation in Graphs and the Global Markov Property . . . 358
14.4.2 The Global Granger-Causal Markov Property . . . . . . . . 360
14.4.3 Markov Properties for Bivariate Path Diagrams . . . . . . . 363
14.4.4 Comparison of Bivariate and Multivariate Granger Causality 364
14.5 Statistical Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . 366
14.5.1 Inference in the Time Domain . . . . . . . . . . . . . . . . . 366
14.5.2 Inference in the Frequency Domain . . . . . . . . . . . . . . 367
14.5.3 Graphical Modeling . . . . . . . . . . . . . . . . . . . . . . . 368
14.6 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 369
14.6.1 Frequency-Domain Analysis of Multivariate Time Series . . 370
14.6.2 Identification of Tremor-Related Pathways . . . . . . . . . . 375
14.6.3 Causal Inference . . . . . . . . . . . . . . . . . . . . . . . . . 377
14.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 379
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 380
15 Multivariate Signal Analysis by Parametric Models
(K. J. Blinowska and M. Kamiñski) 387
15.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 388
15.2 Parametric Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . 388
15.3 Linear Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 390
15.4 Model Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 391
15.5 Cross Measures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393

Contents 7
15.6 Causal Estimators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 394
15.7 Modeling of Dynamic Processes . . . . . . . . . . . . . . . . . . . . 396
15.8 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 398
15.8.1 Common Source in Three Channels System . . . . . . . . . 398
15.8.2 Activity Sink in Five Channels System . . . . . . . . . . . . 398
15.8.3 Cascade Flows . . . . . . . . . . . . . . . . . . . . . . . . . . 402
15.8.4 Comparison between DTF and PDC . . . . . . . . . . . . . . 406
15.9 Multivariate Analysis of Experimental Data . . . . . . . . . . . . . . 408
15.9.1 Human Sleep Data . . . . . . . . . . . . . . . . . . . . . . . . 408
15.9.2 Application of a Time-Varying Estimator of Directedness . 414
15.10Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 417
15.11Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . 420
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 420
16 Computer Intensive Testing for the Influence Between Time Series
(Luiz A. Baccalá, Daniel Y. Takahashi, and Koichi Sameshima) 425
16.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 425
16.2 Basic Resampling Concepts . . . . . . . . . . . . . . . . . . . . . . . 428
16.3 Time Series Resampling . . . . . . . . . . . . . . . . . . . . . . . . . 429
16.3.1 Residue Resampling . . . . . . . . . . . . . . . . . . . . . . . 431
16.3.2 Phase Resampling . . . . . . . . . . . . . . . . . . . . . . . . 432
16.3.3 Other Resampling Methods . . . . . . . . . . . . . . . . . . . 434
16.4 Numerical Examples and Applications . . . . . . . . . . . . . . . . 434
16.4.1 Simulated Data . . . . . . . . . . . . . . . . . . . . . . . . . . 434
16.4.2 Real Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 440
16.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444
16.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 447
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 447
17 Granger Causality: Basic Theory and Application to Neuroscience
(Mingzhou Ding, Yonghong Chen, and Steven L. Bressler) 451
17.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 451
17.2 Bivariate Time Series and Pairwise Granger Causality . . . . . . . . 452
17.2.1 Time Domain Formulation . . . . . . . . . . . . . . . . . . . 452
17.2.2 Frequency Domain Formulation . . . . . . . . . . . . . . . . 454
17.3 Trivariate Time Series and Conditional Granger Causality . . . . . 457
17.3.1 Time Domain Formulation . . . . . . . . . . . . . . . . . . . 458
17.3.2 Frequency Domain Formulation . . . . . . . . . . . . . . . . 459
17.4 Estimation of Autoregressive Models . . . . . . . . . . . . . . . . . 461
17.5 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . 463
17.5.1 Example 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 463
17.5.2 Example 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 465
17.5.3 Example 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 466
17.6 Analysis of a Beta Oscillation Network in Sensorimotor Cortex . . 468

8 Contents
17.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473
18 Granger Causality on Spatial Manifolds: Applications to Neuroimaging
(P. A. Valdés-Sosa, J.M. Bornot-Sánchez, M. Vega-Hernández, L. Melie-García,
A. Lage-Castellanos, and E. Canales-Rodríguez) 475
18.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 476
18.2 The Continuous Spatial Multivariate Autoregressive Model and
Its Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 478
18.3 Testing for Spatial Granger Causality . . . . . . . . . . . . . . . . . 480
18.4 Dimension Reduction Approaches to sMAR Models . . . . . . . . . 482
18.4.1 ROI-Based Causality Analysis . . . . . . . . . . . . . . . . . 482
18.4.2 Latent Variable-Based Causality Analysis . . . . . . . . . . . 483
18.5 Penalized sMAR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 485
18.5.1 General Model . . . . . . . . . . . . . . . . . . . . . . . . . . 485
18.5.2 Achieving Sparsity Via Variable Selection . . . . . . . . . . . 488
18.5.3 Achieving Spatial Smoothness . . . . . . . . . . . . . . . . . 490
18.5.4 Achieving Sparseness and Smoothness . . . . . . . . . . . . 491
18.6 Estimation via the MM Algorithm . . . . . . . . . . . . . . . . . . . 492
18.7 Evaluation of Simulated Data . . . . . . . . . . . . . . . . . . . . . . 494
18.8 Influence Fields for Real Data . . . . . . . . . . . . . . . . . . . . . . 496
18.9 Possible Extensions and Conclusions . . . . . . . . . . . . . . . . . . 498
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 499
Index 505

10 Synchronization Analysis and Recurrence
in Complex Systems
M. C. Romano, M. Thiel, Jürgen Kurths, M. Rolfs, R. Engbert, and R. Kliegl
Author: Please check
your address care-
fully. Please provide
complete address
(full Christian and
family name, insti-
tutional affiliation,
correct address incl.
zip-code, and e-mail)
to be inserted into the
List of Contributors
at the beginning of
the book.
We discuss an approach to detect and quantify phase synchronization in the
case of coupled non-phase-coherent oscillators, which is based on the recurrence
properties of the underlying system. First, we present an index which detects
phase synchronization without computing the phase directly. We show that this
index is also appropriate for non-phase-coherent systems, i.e., systems with a
rather broad power spectrum. Furthermore, we illustrate the applicability of this
index for time series strongly contaminated by noise.
Second, we present an algorithm, which is also based on recurrence to gen-
erate surrogates to test for phase synchronization. The generated surrogates cor-
respond to independent copies of the underlying system. Hence, computing a
phase synchronization index between one observed oscillator and the surrogate
of the second oscillator, we can test for phase synchronization.
Finally, we apply the recurrence-based index, as well as the recurrence-based
surrogates to fixational eye movements and find strong indications that both the
left and right fixational eye movements are synchronized.
10.1 Introduction
The study of synchronization goes back to the seventeenth century and begins
with the analysis of synchronization of nonlinear periodic systems. The synchro-
nization phenomenon was probably discovered first by Huygens in 1673, who
noticed that two pendulum clocks that hang on the same beam can synchronize.
This discovery can be considered as the beginning of Nonlinear Science. The syn-
chronization of the flashing of fireflies, the peculiarities of adjacent organ pipes
which can almost annihilate each other or speak in unison, or the synchroniza-
tion of diodes are other well known examples.
However, the research of synchronization in complex systems did not be-
gin until the end of the eighties. It has been studied extensively during the last
years [1–4], as this phenomenon has found numerous applications in natural
(cardiorespiration, Parkinson patients, ecology, El Niño-Monsoon, etc.) [5–10]
and engineering (lasers, plasma, tubes, etc.) systems [11–13]. Two systems are
Handbook of Time Series Analysis. Björn Schelter, Matthias Winterhalder, Jens Timmer
Copyright © 2006 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim
ISBN: ???????????????

244 10 Synchronization Analysis and Recurrence in Complex Systems
said to be phase synchronized when their respective frequencies and phases are
locked. Note that synchronization is a process (of adapting rhythms) and not a
state. Till now phase synchronization (PS) of chaotic systems has been mainly
observed for attractors with rather coherent phase dynamics. These attractors
have a relatively simple topology of oscillations and a well-pronounced peak in
the power spectrum, which allows to introduce the phase and the characteristic
frequency of motions, Eq. (10.2). However, some difficulties appear when dealing
with non-coherent attractors characterized by a rather broad band power spec-
tra. Then it might not be straightforward to define a phase of the oscillations,
and in general no single characteristic time scale exists. In contrast to phase co-
herent attractors, it is quite unclear whether some phase synchronized state can
be achieved (Fig. 10.1).
To treat this problem, we propose a method based on another basic property
of complex chaotic systems: recurrences in phase space. The concept of recur-
rence in dynamical systems goes back to Poincaré [14], when he proved that after
a sufficiently long time interval, the trajectory of an isolated mechanical system
will return arbitrarily close to each former point of its route. We will show that
the concept of recurrence allows to detect indirectly synchronization and works
even in the case of noisy non-phase-coherent oscillators. Instead of defining di-
rectly the phase, we consider the coincidence of certain recurrence structures of
both coupled subsystems. By means of this comparison we are able to detect syn-
chronization in complex systems.
Another important problem in the synchronization analysis is that even
though the synchronization measures may be normalized, experimental time se-
ries often yield values which are not at the borders of the interval and hence
are difficult to interpret. This problem can be overcome if the coupling strength
between the two systems can be varied systematically and a rather large change
in the measure can be observed, i.e., we have a so called active experiment [1–4].
However, there are other kind of experiments (passive ones), in which it is not
possible to change the coupling strength systematically, e.g., the synchronization
of the heart beats of a mother with her fetus [15]. In some cases, this problem has
been tackled by interchanging the pairs of oscillators [15], for example the EEGs
of other pregnant women were used as “natural surrogates.” These surrogates
are independent and hence not in PS with the original system. Hence, if the syn-
chronization index obtained with the original data is not significantly higher than
the index obtained with the natural surrogates, there is no sufficient evidence to
claim synchronization. But even this rather innovative approach has some draw-
backs. The natural variability and also the frequency of the heart beats of the
surrogate mothers are usually lightly different from the ones of the real mother.
Furthermore, the data acquisition can be expensive and at least in some cases
problematic or even impossible (e.g., geophysical time series). In these cases it
would be convenient to perform a hypothesis test based on surrogates generated
by a mathematical algorithm.

10.2 Phase Synchronization by Means of Recurrences 245
Therefore, we present a technique for the generation of surrogates, which is
based on the recurrences of a system. These surrogates mimic the dynamical
behavior of the system. Then, computing the synchronization index between one
subsystem of the original system and the other subsystem of the surrogate, and
comparing it with the synchronization index obtained for the original system, we
can test for PS.
In Section 10.2, we introduce the concept of recurrence, as well as the synchro-
nization index based on the recurrence properties of the system. In Section 10.3
we show how to detect another kind of synchronization, namely generalized
synchronization (GS) by means of recurrences and in Section 10.4 we show that
the recurrence-based indices indicate the transition to PS and GS in accordance
with other known theoretical methods. In Section 10.5 we present the twin sur-
rogates technique and apply it to test for synchronization in the paradigmatic
two coupled Rössler systems. In Section 10.6 we show an application of the
recurrence-based index and surrogates to measured physiological data, namely
fixational eye movements.
10.2 Phase Synchronization by Means of Recurrences
First, we exemplify the problem of defining the phase in systems with rather
broad power spectrum by the paradigmatic system of two coupled nonidentical
Rössler oscillators
x˙1,2 = −ω1,2y1,2 − z1,2 ,
y˙1,2 = ω1,2x1,2 + ay1,2 + µ(y2,1 − y1,2) ,
z˙1,2 = 0.1 + z1,2(x1,2 − 8.5) ,
(10.1)
where µ is the coupling strength and ω1,2 determine the mean intrinsic fre-
quency of the (uncoupled) oscillators in the case of phase coherent attractors. In
our simulations we take ω1 = 0.98 and ω2 = 1.02. The parameter a ∈ [0.15 : 0.3]
governs the topology of the chaotic attractor. When a is below a critical value ac
(ac ≈ 0.186 for ω1 = 0.98 and ac ≈ 0.195 for ω2 = 1.02), the chaotic trajectories
always cycle around the unstable fixed point (x0, y0) ≈ (0, 0) in the (x, y) sub-
space, i.e., max(y) > y0 (Fig. 10.1(a)). In this case, simply the rotation angle
φ = arctan
y
x
(10.2)
can be defined as the phase, which increases almost uniformly. The oscillator has
a coherent phase dynamics, i.e., the diffusion of the phase dynamics is very low
(10−5 to 10−4). In this case, other phase definitions, e.g., based on the Hilbert
transform or on the Poincaré section, yield equivalent results [1–4]. However,
beyond the critical value ac, the trajectories no longer always completely cycle
around (x0, y0), and some max(y) < y0 occur, which are associated with faster
returns of the orbits; the attractor becomes a funnel one. Such earlier returns in

246 10 Synchronization Analysis and Recurrence in Complex Systems
Fig. 10.1: (a,c) Segment of the x1-component of the trajectory of the Rössler sys-
tems, Eq. (10.1). (b,d) periodogram of the x-component of the trajectory. (e,g) pro-
jection of the attractor onto the (x, y) plane. (g,h) projection onto the (x˙, y˙) plane.
(a,b,e,f) computed for a = 0.16 and (c,d,g,h) computed for a = 0.2925.
the funnel attractor happen more frequently with increasing a (Fig. 10.1(b)). It is
clear that for the funnel attractors, usual (and rather simple) definitions of phase,
such as Eq. (10.2), are no longer applicable [1–4].
Rosenblum et al. have proposed in [16] to use an ensemble of phase coherent
oscillators which is driven by the non-phase-coherent oscillator in order to esti-
mate the frequency of the last, and hence detect PS in such kind of systems. How-
ever, depending on the component one uses to couple the non-phase-coherent
oscillator to the coherent ones, the result of the obtained frequency can be differ-
ent.
Furthermore, Osipov et al. [17] have proposed another approach which is
based on the general idea of the curvature of an arbitrary curve [18]. For any two-
dimensional curve r1 = (u, v) the angle velocity at each point is ν =
ds
dt/R, where
ds/dt =
√
u˙2 + v˙2 is the velocity along the curve and R = (u˙2 + v˙2)3/2/[v˙u¨ −
v¨u˙] is the curvature. If R > 0 at each point, then ν = dφdt =
v˙u¨−v¨u˙
u˙2+v˙2 is always
positive and therefore the variable φ defined as φ =
∫
νdt = arctan v˙u˙ , is a
monotonically growing angle function of time and can be considered as a phase

10.2 Phase Synchronization by Means of Recurrences 247
of the oscillations. Geometrically it means that the projection r2 = (u˙, v˙) is a
curve cycling monotonically around a certain point.
These definitions of φ and ν hold in general for any dynamical system if
the projection of the phase trajectory onto some plane is a curve with a positive
curvature. This approach is applicable to a large variety of chaotic oscillators,
such as the Lorenz system [19], the Chua circuit [20] or the model of an ideal
four-level laser with periodic pump modulation [21].
This is clear for phase-coherent as well as funnel attractors in the Rössler
oscillator. Here projections of chaotic trajectories on the plane (x˙, y˙) always rotate
around the origin (Fig. 10.1(c) and (d)) and the phase can be defined as
φ = arctan
y˙
x˙
. (10.3)
We have to note that for funnel-like chaotic attractors the coupling may change
their topology. As a consequence the strong cyclic structure of orbits projection
in the (x˙, y˙)-plane may be destroyed and the phase measurement by Eq. (10.3)
fails occasionally for intermediate values of coupling. But for small coupling and
for coupling near the transition to PS, the phase is well-defined by Eq. (10.3) [22].
We consider two criteria to detect the existence of PS: Locking of the mean
frequencies Ω1 = 〈ν1〉 = Ω2 = 〈ν2〉, and locking of the phase |φ2(t) − φ1(t)|

const (we restrict here to 1 : 1 synchronization). Applying the new definition of
the phase Eq. (10.3) to the system defined by Eq. (10.1) for a = 0.2925 (strongly
noncoherent) and µ = 0.179, one obtains the phase difference represented in
Fig. 10.2.
We find two large plateaus in the evolution of the difference of the phases
with time, i.e., we detect PS, but we also find a phase slip associated to a dif-
ferent number of oscillations in the two oscillators in the represented period of
time. This means, we observe the rare occurrence of phase slip. It is interesting to
note that in this system PS occurs after one of the positive Lyapunov exponents
passes to negative values, i.e., it is also a transition to generalized chaotic syn-
chronization (GS).
Although this approach works well in non-phase-coherent model systems, we
have to consider that one is often confronted with the computation of the phase
in experimental time series, which are usually corrupted by noise. In this case,
some difficulties may appear when computing the phase by Eq. (10.3), because
derivatives are involved in its definition.
Hence, we propose a different approach based on recurrences in phase space
to detect PS indirectly. We define a recurrence of the trajectory of a dynamical
system {x(i)}Ni=1 in the following way: We say that the trajectory has returned at
time t = jδt to the former point in phase space visited at t = iδt if
R(ε)i,j = Θ(ε − ‖x(i) − x(j)‖) = 1 , (10.4)
where ε is a pre-defined threshold, Θ(.) is the Heaviside function and δt is the
sampling rate. A “1” in the matrix at i, j means that x(i) and x(i) are neighbor-

248 10 Synchronization Analysis and Recurrence in Complex Systems
27100 27150 27200 27250 27300 27350
time
−20
−10
0
10
20
dy
1/d
t,
dy
2/d
t
−2
2
6
10
φ 2−
φ 1
(a)
(b)
Fig. 10.2: (a) Time evolution of phase difference of the system of Eq. (10.1). (b)
Variables y˙1,2 in system (10.1) for a = 0.2925 and µ = 0.179. Solid and dotted
lines correspond to the first and the second oscillator, respectively. In the time
interval between dashed lines the first oscillator produces four rotations in the
(x˙1, y˙1)-plane around the origin, but the second one generates only three rotations,
which leads to a phase slip in (a).
ing, a “0” that they are not. The black and white representation of this binary
matrix is called recurrence plot (RP). This method has been intensively studied
in the last years: Different measures of complexity have been proposed based on
the structures obtained in the RP and have found numerous applications in, e.g.,
physiology and earth science [23–27]. Furthermore, it has been even shown that
some dynamical invariants can be estimated by means of the recurrence struc-
tures [28].
Based on this definition of recurrence, we want to tackle the problem of per-
forming a synchronization analysis in the case of non-phase-coherent systems.
We avoid the direct definition of the phase and instead use the recurrence prop-
erties of the systems in the following way: The probability P(ε)(τ) that the system
returns to the neighborhood of a former point x(i) of the trajectory1 after τ time
steps can be estimated as follows:
Pˆ(ε)(τ) =
1
N − τ
N−τ
∑
i=1
Θ(ε − ‖x(i) − x(i + τ)‖) = 1
N − τ
N−τ
∑
i=1
R(ε)i,i+τ . (10.5)
This function can be regarded as a generalized autocorrelation function, as it
also describes higher order correlations between the points of the trajectory in
1 The neighborhood is defined as a box of size ε centered at x(i), as we use the maximum norm.

10.2 Phase Synchronization by Means of Recurrences 249
dependence on the time delay τ. A further advantage with respect to the linear
autocorrelation function is that Pˆ(ε)(τ) is defined for a trajectory in phase space
and not only for a single observable of the system’s trajectory. Further, we have
recently shown that it is possible to reconstruct the attractor by only considering
the recurrences of single components of the system [29] and it is also possible to
estimate dynamical invariants of the system (e.g., entropies and dimensions) by
means of recurrences in phase space [28]. Hence, the recurrences of the system
in phase space contain information about higher order dependencies within the
components of the system.
For a periodic system with period length T in a two-dimensional phase space,
it can be easily shown that
P(τ) = lim
ε→0
Pˆ(ε)(τ) =
{
1 τ = T
0 otherwise .
For coherent chaotic oscillators, such as Eq. (10.1) for a = 0.16, Pˆ(ε)(τ) has well-
expressed local maxima at multiples of the mean period, but the probability of
recurrence after one or more rotations around the fixed point is less than one
(Fig. 10.5).
Analyzing the probability of recurrence, it is possible to detect PS for non-
phase-coherent oscillators, too. This approach is based on the following idea:
Originally, a phase φ is assigned to a periodic trajectory x in phase space, by
projecting the trajectory onto a plane and choosing an origin, around which the
trajectory oscillates all the time. Then an increment of 2π is assigned to φ, when
the point of the trajectory has returned to its starting position, i.e., when ‖x(t +
T) − x(t)‖ = 0. Analogously to the case of a periodic system, we can refer an
increment of 2π to φ to a complex nonperiodic trajectory x(t), when ‖x(t + T) −
x(t)‖ ∼ 0, or equivalently when ‖x(t + T) − x(t)‖ < ε, where ε is a predefined
threshold. That means, a recurrence R(ε)t,t+τ = 1 can be interpreted as an increment
of 2π of the phase in the time interval τ2.
Pˆ(ε)(τ) can be viewed as a statistical measure on how often φ in the original
phase space has increased by 2π or multiples of 2π within the time interval τ. If
two systems are in PS, in the mean, the phases of both systems increase by k · 2π,
with k a natural number, within the same time interval τ. Hence, looking at the
coincidence of the positions of the maxima of Pˆ(ε)(τ) for both systems, we can
quantitatively identify PS (from now on, we omit (ε) and ·ˆ in Pˆ(ε)(τ) to simplify
the notation). The proposed algorithm then consists of two steps:
• Compute P1,2(τ) of both systems based on Eq. (10.5).
• Compute the cross-correlation coefficient between P1(τ) and P2(τ) (correlation
between probabilities of recurrence, CPR)
CPR1,2 =
〈P¯1(τ)P¯2(τ)〉τ
σ1σ2
, (10.6)
2 This can be considered as an alternative definition of the phase to Eqs. (10.2) and (10.3).

250 10 Synchronization Analysis and Recurrence in Complex Systems
Fig. 10.3: P(τ) for a periodically driven Rössler (Eqs. (10.7)) in PS (a) and in non-PS
(b). Solid line: P(τ) of the driven Rössler, dashed line: P(τ) of the periodic forcing.
where P¯1,2 means that the mean value has been subtracted and σ1 and σ2 are
the standard deviations of P1(τ) and P2(τ), respectively.
If both systems are in PS, the probability of recurrence is maximal simultane-
ously and CPR1,2 ∼ 1. In contrast, if the systems are not in PS, the maxima of the
probability of recurrence do not occur jointly and expect low values of CPR1,2.
10.2.1 Examples of Application
In this section we exemplify the application of CPR to detect PS for four prototyp-
ical chaotic systems. The number of data points used for the analysis presented
here is 5000.
1. We start with the periodically driven Rössler system [1–4]
x˙ = −y − z + µ cos(ωt)
y˙ = x + 0.15y
z˙ = 0.4 + z(x − 8.5) .
(10.7)
For the frequency ω = 1.04 and the coupling strength µ = 0.16, the periodic
forcing locks the frequency of the Rössler system. This can be clearly seen in
Fig. 10.3(a). The position of the maxima coincide. The value of the recurrence-
based PS index (Eq. (10.6)) is CPR = 0.862.
For the parameters ω = 1.1 and µ = 0.16, the periodic forcing does not syn-
chronize the Rössler system. Hence, the joint probability of recurrence is very
low, which is reflected in the drift of the peaks of the corresponding P(τ)
(Fig. 10.3(b)). In this case, CPR = −0.00241.

10.2 Phase Synchronization by Means of Recurrences 251
Fig. 10.4: P(τ) for a periodically driven Lorenz in PS (a) and in non-PS (b). Solid
line: P(τ) of the driven Lorenz, dashed line: P(τ) of the periodic forcing.
2. We continue our considerations with the periodically driven Lorenz system
for the standard parameters
x˙ = 10(y − x)
y˙ = 28x − y − xz
z˙ = −8/3z + xy + µ cos(ωt) .
(10.8)
In Fig. 10.4(a) the probabilities of recurrence P(τ) in the PS case (µ = 10,
ω = 8.35) are represented. We see that the position of the local maxima of the
Lorenz oscillator coincide with the ones of the periodic forcing. However, the
local maxima are not as high as in the case of the Rössler system, and they
are broader. This reflects the effective noise which is intrinsic in the Lorenz
system [1–4]. Therefore, the phase synchronization is not perfect: An exact
frequency locking between the periodic forcing and the driven Lorenz cannot
be observed [30]. In this case, we obtain CPR = 0.667. In the non-PS case
(µ = 10, ω = 7.5), we obtain CPR = 0.147 (Fig. 10.4(b)).
3. Next, we consider the case of two mutually coupled Rössler systems in the
phase coherent regime, more precisely we analyze Eqs. (10.1) with a = 0.16.
According to [17], for ω1 = 0.98,ω2 = 1.02 and µ = 0.05 both systems are in
PS. We observe that the local maxima of P1 and P2 occur at τ = n · T , where
T is the mean period of both Rössler systems (Fig. 10.5(a)). The heights of
the local maxima are in general different for both systems if they are only in
PS and not in, e.g., complete synchronization or generalized synchronization.
But the positions of the local maxima of P(τ) coincide. In this case, we obtain
CPR = 0.998.

252 10 Synchronization Analysis and Recurrence in Complex Systems
Fig. 10.5: P(τ) for two mutually coupled Rössler systems (Eqs. (10.1)) in phase
coherent regime (a = 0.16) for µ = 0.05 (a) and for µ = 0.02 (b).
At a coupling strength of µ = 0.02 the systems are not in PS and the positions
of the maxima of P(τ) do not coincide anymore (Fig. 10.5(b)), clearly indicating
that the frequencies are not locked. In this case, we obtain CPR = 0.115.
4. As a last example with simulated data, we analyze the challenging case of two
mutually coupled Rössler systems in the funnel regime. Therefore, we study
Eqs. (10.1) with a = 0.2925, ω1 = 0.98, and ω2 = 1.02. We analyze two dif-
ferent coupling strengths: µ = 0.2 and µ = 0.05. We observe that the structure
of P(τ) in the funnel regime (Fig. 10.6) is rather different from the one in the
phase coherent Rössler system (Fig. 10.5). The peaks in P(τ) are not as well
pronounced as in the coherent regime, reflecting the different time scales that
play a crucial role and the broad-band power spectrum of this system. How-
ever, we notice that for µ = 0.2 the locations of the local maxima coincide for
both oscillators (Fig. 10.6(a)), indicating PS, whereas for µ = 0.05 the positions
of the local maxima do not coincide anymore (Fig. 10.6(b)), indicating non-PS.
These results are in accordance with [17].
In the PS case, we obtain CPR = 0.988, and in the non-PS case, CPR = 0.145.
Note that the position of the first peak in Fig. 10.6(b) coincides, although the
oscillators are not in PS. This is due to the small frequency mismatch (|ω1 −
ω2| = 0.04). However, by means of the index CPR we can distinguish rather
well between both regimes.

10.2 Phase Synchronization by Means of Recurrences 253
Fig. 10.6: P(τ) for two mutually coupled Rössler systems (Eqs. (10.1)) in funnel
regime (a = 0.2925) for µ = 0.2 (a) and for µ = 0.05 (b). Bold line: P1(τ), solid line:
P2(τ).
Fig. 10.7: First component x1 of Eqs. (10.1) with 80% independent Gaussian noise
(for µ = 0.05). From the figure it is clearly recognizable that it is difficult to com-
pute the phase by means of, e.g., the Hilbert transformation.
10.2.2 Influence of Noise
Measurement errors are omnipresent in experimental time series. Hence, it is
necessary to analyze the influence of noise on CPR (correlation of probability of
recurrence).

254 10 Synchronization Analysis and Recurrence in Complex Systems
Fig. 10.8: Probabilities of recurrence for two coupled Rössler systems (Eqs. (10.1))
in PS (µ = 0.05) without noise (a) and with 80% Gaussian observational noise (b).
Bold line: subsystem 1, solid line: subsystem 2. Note that the position of the peaks
of P1(τ) and P2(τ) coincide in both cases, and hence the solid line is hidden by the
bold one.
First, we treat additive or observational white noise. We use Eqs. (10.1) with
two different coupling strengths, so that we can compute the deviation which is
caused by noise in the nonsynchronized and in the synchronized case.
We add independent Gaussian noise with standard deviation σnoise = ασj
to each coordinate j of the system, where σj is the standard deviation of the
component j and α is the noise level. In Fig. 10.7 the “corrupted” x-component
of the first Rössler subsystem x˜1(t) = x1(t) + ασ1η(t) is represented. Herein η(t)
is a realization of Gaussian noise and α = 0.8. From Fig. 10.7 it is obvious that it
is difficult to compute the phase by means of, e.g., the Hilbert transformation for
such a high noise level without filtering.
The choice of ε for the computation of P1(τ) and P2(τ) in the presence of noise
is performed automatically by fixing the recurrence rate RR, i.e., the percentage
of recurrence points in the recurrence matrix, Eq. (10.4). The results presented
below were computed for RR = 0.1, but the results are rather independent of the
choice of RR. However, RR should not be chosen too small if the level of noise is
very high [23–27].
In order to compute CPR for the noisy oscillators, we calculate first the prob-
abilities of recurrence P1(τ) and P2(τ) for coupling strengths µ = 0.05 (PS,
Fig. 10.8) and µ = 0.02 (non-PS, Fig. 10.9).
We note that the peaks in P1(τ) and P2(τ) become lower and broader (Figs. 10.8(b)
and 10.9(b)) compared with the noise free case (Figs. 10.8(a) and 10.9(a)), which
is expected. However, despite of the large level of noise, the positions of the local

10.2 Phase Synchronization by Means of Recurrences 255
Fig. 10.9: Probabilities of recurrence for two coupled Rössler systems (Eqs. (10.1))
in non-PS (µ = 0.02) without noise (a) and with 80% Gaussian observational
noise (b). Bold line: subsystem 1, solid line: subsystem 2.
maxima coincide in the PS case, and they drift away in the non-PS case. This
a convenient result, because we can still decide whether the oscillators are syn-
chronized in a statistical sense or not. This is reflected in the obtained values for
the CPR index: at a noise level of 80% noise, in the PS case the obtained value
for CPR is exactly the same with and without noise, and in the non-PS case it is
nearly the same (see Table 10.1). This shows that the index CPR for PS is very
robust against observational noise.
Now, we analyze the influence of colored noise on the index CPR. We add a
realization of colored noise with a very high noise amplitude to each component
of the first system and another realization of colored noise with a smaller noise
amplitude to each component of the second system (see Fig. 10.10(a) and (b) and
the corresponding caption). Other methods fail determining the phase in this
case, as for example the one presented in [17], because it requires the computation
of the derivative of the time series, and due to the large level of noise, this is not
possible. But by means of P(τ) we can distinguish PS from non-PS even in this
case (Fig. 10.10(c) and (d)): We obtain CPR = 0.0276 for the non-PS case and
CPR = 0.530 for the PS case.
Tab. 10.1: Index CPR for PS calculated for two coupled Rössler systems (10.1) with
observational noise and without noise, for comparison.
µ CPR (80% noise) CPR (0% noise)
0.02 (non-PS) 0.149 0.115
0.05 (PS) 0.998 0.998

256 10 Synchronization Analysis and Recurrence in Complex Systems
Fig. 10.10: (a,b) Segments of the x-components of the trajectories of two mutually
coupled Rössler systems in phase coherent regime (a = 0.16) strongly contami-
nated by colored noise. A realization of rt+1 = 0.99rt+10ηt and st+1 = 0.982st+ξt
were respectively added to each component of the Rössler systems. (a) non-PS
(µ = 0.02). (b) PS (µ = 0.05). (c) P(τ) for the two noisy Rössler for µ = 0.02
(non-PS), (d) P(τ) for the two noisy Rössler for µ = 0.05 (PS). Solid line: system 1,
dashed line: system 2.
10.3 Generalized Synchronization and Recurrence
In this section we treat the issue of synchronization of coupled systems which
are essentially different. This problem has been addressed first in [31, 32]. In this
case, there is in general no trivial manifold in the phase space which attracts
the systems’ trajectories. It has been shown that these systems can synchronize
in a more general way, namely y = ψ(x), where ψ is a transformation which
maps asymptotically the trajectories of x into the ones of the attractor y. This
kind of synchronization is called generalized synchronization (GS). The proper-
ties of the function ψ depend on the features of the systems x and y, as well as
on the attraction properties of the synchronization manifold y = ψ(x) [33]. GS
has been demonstrated in laboratory experiments for electronic circuits and laser
systems [34–38] and has found applications for the the design of communication
devices [39–43] and model verification and parameter estimations from time se-
ries [44, 45].
Some statistical measures have been introduced for the detection of GS, such
as the method of mutual false nearest neighbors [31, 32] or variations of the
method proposed and analyzed in [46–48], which are based on the squared mean
distance and conditional distance between mutual nearest neighbors. Some other
methods are based on the mutual predictability to detect dynamical interdepen-
dence [49, 50]. There, the nearest neighbors of each subsystem are computed
separately in the respective (sub)state space.

10.3 Generalized Synchronization and Recurrence 257
In this section we present a criterion for the detection of GS, which exploits the
relationship between the geometric connection of both systems and their recur-
rences. The connection between recurrences and GS is even more straightforward
than the one between recurrences and PS. One can see that the concept of GS is
linked to the one of recurrence, considering the fact that when x(t) and y(t) are
in GS, two close states in the phase space of x correspond to two close states
in the space of y [31, 32]. Hence, the “neighborhood identity” in phase space
is preserved, i.e., they are topologically equivalent. Since the recurrence matrix
(Eq. (10.4)) is nothing else but a record of the neighborhood of each point of the
trajectory, one can conclude that two systems are in GS if their respective RPs are
almost identical. Note that it is possible, under some conditions, to reconstruct
the rank order of the time series considering only the information contained in
the RP [29]. Therefore, we can use the recurrence properties to detect and quan-
tify GS.
However, in practice we note that the recurrence matrices of two systems in
GS are very similar, but not identical. Several reasons can be given to explain
this observation: The finite ε-threshold, computational roundoff errors, measure-
ments inaccuracies, etc. Hence, we construct an index that quantifies the degree
of similarity between the respective recurrences of both systems. It compares the
recurrences of each point of the first system with the local recurrences of the
second system. This index has the advantage that it distinguishes rather well be-
tween non-PS, PS, and GS.
This index is based on the average probability of joint recurrence over time,
given by
RRx,y =
1
N2
N
∑
i,j=1
Θ(εx − ‖xi − xj‖)Θ(εy − ‖yi − yj‖) . (10.9)

258 10 Synchronization Analysis and Recurrence in Complex Systems
If both systems x and y are independent from each other, then the average
probability of a joint recurrence3 is given by RRx,y = RRxRRy. If the oscillators
are on the other hand in GS, we expect an approximate identity of their respective
recurrences, and hence RRx,y = RRx = RRy [31, 32].
For the computation of the recurrence matrix in the case of essentially dif-
ferent systems that undergo GS, it is more appropriate to use a fixed number
of nearest neighbors for each column in the matrix, following the idea pre-
sented in [46–48], than using a fixed threshold. This means that the threshold
is different for each column in the RP, but subjected to the following condition
∑N
j=1 Θ(ε
i −‖xi −xj‖) = A ∀i, where A is the fixed number of nearest neighbors.
We can automatically fix the RR by means of RR = AN/N2 = A/N, and using
the same A for each subsystem x and y, RRx = RRy = RR.
Hence, the coefficient S = RR
x,y
RR is an index for GS that varies from RR to 1:
It is approximately RR for independent systems, and it is close to 1 for systems
in GS. However, with the index S we would not detect lag synchronization (LS)
(y(t+τ) = x(t)). Since LS can be considered as a special case of GS [52], it would
be desirable to have an index that also detects LS. For this reason, we include a
time lag τ in the similarity and introduce the following quotient:
S(τ) =
1/N2
∑N
i,j Θ(ε
i
x − lVertxi − xj‖)Θ(εiy − ‖yi+τ − yj+τ‖)
RR
, (10.10)
where the thresholds εix and ε
i
y fullfil the following conditions:
∑N
j=1 Θ(ε
i
x−‖xi−
xj‖) = A and
∑N
j=1 Θ(ε
i
y − ‖yi − yj‖) = A ∀i. Then, we choose the maximum
value of S(τ) and normalize
JPR = max
τ
S(τ) − RR
1 − RR
. (10.11)
We denote this index by JPR because it is based on the average joint probability
of recurrence. Since S(τ) varies between RR and 1, JPR ranges from 0 to 1. The
value of RR is a free parameter and its choice depends on the case under study.
We consider rather low values of RR, e.g., 1% or 2% as appropriate.
10.3.1 Examples of Application
In this section we show two examples of chaotic systems that undergo GS and
compute for them the recurrence-based index JPR (Eq. (10.11)).
1. First we consider a Lorenz system driven by a Rössler system. The equations
of the driving system are
x˙1 = 2 + x1(x2 − 4)
x˙2 = −x1 + x3
x˙3 = x2 + 0.45x3 ,
(10.12)
3 Note that the average probability of a joint recurrence is the recurrence rate of the joint recurrence
plot (JRP) [51].

10.3 Generalized Synchronization and Recurrence 259
Fig. 10.11: Projection of the Rössler driving system (a), the driven Lorenz system (b)
and the diagram x2 versus y2 of Eqs. (10.12) and (10.13) (c).
and following are the equations of the driven system:
y˙1 = −σ(y1 − y2)
y˙2 = ru(t) − y2 − u(t)y3
y˙3 = u(t)y2 − by3 ,
(10.13)
where u(t) = x1(t) + x2(t) + x3(t) and the parameters were chosen as follows:
σ = 10, r = 28, and b = 2.666. In [53] it was shown that the systems given by
Eqs. (10.12) and (10.13) are in GS, since the driven Lorenz system is asymptot-
ically stable.
To illustrate that they are completely different systems and that they are not
in LS or even complete synchronized, Fig. 10.11 shows the projections of the
system (Eqs. (10.12)) (a), of the system (Eqs. (10.13)) (b) and the x2 versus y2
diagram (c).
When dealing with experimental time series, usually only one observable of
the system is available. Hence, we perform the analysis with just one com-
ponent of each system to illustrate the applicability of the proposed method
(we use 10 000 data points with a sampling time interval of 0.02 s). In this ex-
ample, we take x3 and y3 as observables, respectively. Then, we reconstruct
the phase space vectors using delay coordinates [54]. For the subsystem x we
obtain the following embedding parameters [55]: delay time τ = 5 and embed-
ding dimension m = 3. For the subsystem y we find: τ = 5 and m = 7. The
corresponding RPs and JRP are represented in Fig. 10.12.
We see that despite of the essential difference between both subsystems, their
RPs are very similar (Fig. 10.12(a) and (b)). Therefore, the structures are re-
flected also in the JRP and consequently, its recurrence rate is rather high. In
this case, with the choice RR = 0.02 we obtain JPR = 0.605 (the value of JPR is
similar for other choices of RR).
In order to illustrate the second case, where both subsystems are independent
(Fig. 10.13), we compute the RP of the Rössler system (Eqs. (10.12)) and of

260 10 Synchronization Analysis and Recurrence in Complex Systems
Fig. 10.12: (a) RP of the Rössler subsystem (Eqs. (10.12)). (b) RP of the driven
Lorenz subsystem (Eqs. (10.13)). (c) JRP of whole system (Eqs. (10.12) and (10.13)).
the independent Lorenz system,4 as well as their JRP (Fig. 10.14). Note that We placed the refer-
ence in a footnote.
Please check.
the mean probability over time for a joint recurrence is very small, as the JRP
has almost no recurrence points. In this case, one obtains JPR = 0.047 using
embedding parameters τ = 5 and m = 3 for both systems, and RR = 0.02.
For σ = 10 and b = 8/3 they display chaotic behavior.
2. Two mutually coupled Rössler systems (Eqs. (10.1)): for the coupling strength
µ = 0.11 both oscillators are in LS, as can be seen from Fig. 10.15.
In this case, the RPs of both subsystems are obviously almost identical, except
for a displacement on τ in the diagonal direction. Computing the index fol-
lowing Eq. (10.11), we obtain the value JPR = 0.988 (JPR in this case is not
exactly 1), because we do not have perfect LS, i.e., x(t + τ)  y(t) [52]). For a
smaller coupling strength µ = 0.02 the oscillators are not in LS anymore. The
obtained value in this is case JPR = 0.014.
4 The Lorenz equations are given by x˙ = −σx+σy, y˙ = −xz+rx−y, z˙ = xy−bz. For σ = 10
and b = 8/3 they display chaotic behavior.

10.4 Transitions to Synchronization 261
Fig. 10.13: Projection of the Rössler system (Eqs. (10.12)) (a), the independent
Lorenz system (see footnote 4) (b) and the diagram x2 versus y2, where x2 is
the second component of the Rössler system and y2 is the second component of
the independent Lorenz system (c).
Fig. 10.14: (a) RP of the Rössler subsystem (Eqs. (10.12)). (b) RP of the independent
Lorenz system (see footnote 4) (c) JRP of whole system.
10.4 Transitions to Synchronization
We have seen in the previous sections that the indices CPR and JPR clearly distin-
guish between oscillators in PS and oscillators which are not in PS, respectively of

262 10 Synchronization Analysis and Recurrence in Complex Systems
Fig. 10.15: Example of lag synchronization: It is clearly seen that x1 (bold line) goes
behind y1 (solid line). It holds: x1(t + τ) = y1(t), with τ = 4.
GS. On the other hand, the synchronization indices should not only distinguish
between synchronized and nonsynchronized regimes, but also clearly indicate
the onset of PS, respectively of GS.
In order to demonstrate that the recurrence-based indices fulfill this condition,
we exemplify their application in the two cases: Two mutually coupled Rössler
systems in a phase coherent regime and in a non-phase-coherent funnel regime
(Eqs. (10.1)) with a = 0.16, respectively a = 0.2925). In both the cases we increase
the coupling strength µ continuously and compute for each value of µ the indices
CPR and JPR.
On the other hand, in the phase coherent case for a not too large but fixed fre-
quency mismatch between both oscillators and increasing coupling strength, the
transitions to PS and LS are reflected by the Lyapunov spectrum [1–4].5 If both os-
cillators are not in PS, there are two zero Lyapunov exponents (λ3 and λ4), which
correspond to the (almost) independent phases. Increasing the coupling strength,
the fourth Lyapunov exponent λ4 becomes negative (Fig. 10.16(c)), indicating the
onset of PS. For higher coupling strengths, the second Lyapunov exponent λ2
crosses zero, which indicates the establishment of a strong correlation between
the amplitudes (Fig. 10.16(c)). This last transition occurs almost simultaneously
with the onset of LS [52]. Therefore, we compute for our two examples also λ2
and λ4 in order to validate the results obtained with CPR and JPR.
In Fig. 10.16 the indices CPR (a) and JPR (b) are represented for increasing
5 For other cases, e.g., for a fixed coupling strength and decreasing frequency mismatch, or for a
large frequency mismatch and increasing coupling strength, the transition to PS is not always
simply reflected in the Lyapunov spectrum [17, 51].

10.4 Transitions to Synchronization 263
Fig. 10.16: CPR index, JPR index and λ2 and λ4 as functions of the coupling
strength µ for two mutually coupled Rössler systems in phase coherent regime
(a,c,e) and in funnel one (b,d,f). The dotted zero line in (e) and (f) is plotted to
guide the eye. Here, we choose ε corresponding to 10% recurrence points in each
RP.
coupling strength µ for the phase coherent case. In (c) λ2 and λ4 are shown in
dependence on µ.
By means of CPR, the transition to PS is detected when CPR becomes of the
order of 1. We see from Fig. 10.16(a) that the transition to PS occurs at approx-
imately µ = 0.037, in accordance with the transition of the fourth Lyapunov
exponent λ4 to negative values. The index JPR shows three plateaus in depen-
dence on the coupling strength (Fig. 10.16(b)), indicating the onset of PS at the
beginning of the second one. On the other hand, JPR clearly indicates the onset
of LS because it becomes nearly one (third plateau) at approximately µ = 0.1
(Fig. 10.16(b)), after the transition from hyperchaoticity to chaoticity, which takes
place at approximately µ = 0.08 (Fig. 10.16(c)). Between µ = 0.08 and µ = 0.1,
the values of JPR have large fluctuations. This reflects the intermittent LS [1–4],
where LS is interrupted by intermittent bursts of no synchronization.
Now we regard the more complex case of two coupled Rössler systems in the
non-phase-coherent funnel regime, where the direct application of the Hilbert
transformation is not possible [17]. In Fig. 10.16 the coefficients CPR and JPR are
represented for this case in dependence on the coupling strength µ. Again, λ2
and λ4 are also shown (Fig. 10.16(f)).

264 10 Synchronization Analysis and Recurrence in Complex Systems
First, note that for µ > 0.02, λ4 has already passed to negative values (Fig. 10.16(f)).
However, CPR is still rather low, indicating that both oscillators are not in PS yet.
CPR does not indicate PS until µ = 0.18 (Fig. 10.16(d)), as found with other tech-
niques [17]. Furthermore, we see from Fig. 10.16(f) that λ2 vanishes at µ ∼ 0.17.
This transition indicates that the amplitudes of both oscillators become highly
correlated. At approximately the same coupling strength, JPR reaches rather high
values, indicating the transition to GS (Fig. 10.16(e)). Then, according to the in-
dex CPR the transition to PS occurs after the onset of GS. This is a general result
that holds for systems with a strong phase diffusion, as reported in [17]. For
highly non-phase-coherent systems, there is more than one characteristic time
scale. Hence, a rather high coupling strength is necessary in order to obtain phase
locking of both oscillators. Hence, PS is not possible without a strong correlation
in the amplitudes. PS for such non-phase-coherent systems has been recently
found and analyzed in electrochemical oscillators [56] and in El Niño-Monsoon
system [57].
Note that the synchronization indices presented in these sections based on
recurrences are applicable to multivariate time series.
10.5 Twin Surrogates to Test for PS
As we have mentioned in Section 10.1, another essential problem in the synchro-
nization analysis of observed time series is the construction of an appropriate
hypothesis test to test for PS. Several approaches in this direction have been pub-
lished [58, 59]. Usually, these are linear surrogates based on randomization of
the Fourier phases [60, 61]. They mimic the individual spectra of the two com-
ponents of the original bivariate series as well as their cross-spectrum, i.e., their
linear properties, but not the higher order moments. In this case, the correspond-
ing null hypothesis is that the putative synchronization in the underlying system
can be explained by a bivariate linear stochastic process. The specificity of this
test is not always satisfactory, because the concept of PS assumes the mutual
adaption of self-sustained oscillators, i.e., nonlinear deterministic systems. On
the other hand, pseudo-periodic surrogates (PPS) have been proposed to test the
null hypothesis that an observed time series is consistent with an uncorrelated
noise-driven periodic orbit [62]. The PPS are in a certain sense closer to the sur-
rogates needed to test for PS as they correspond to trajectories of a deterministic
system with noise, but they are still not appropriate to test for PS, as they are not
able to model chaotic oscillators. Therefore, we present a technique for the gen-
eration of surrogates which are consistent with the null hypothesis of a trajectory
of the same underlying system, but starting at different initial conditions [63].
Hence, they can also be used to test for PS in the case of chaotic oscillators.
The main idea consists in exchanging one original subsystem with one sur-
rogate. Then, if the synchronization index obtained for the original system is not
significantly different from the one computed for the exchanged subsystems, we

10.5 Twin Surrogates to Test for PS 265
Fig. 10.17: This diagram represents the main idea using twin surrogates to test for
PS.
have no sufficient evidence to claim synchronization (see Fig. 10.17). One could
argue that the same can be achieved using different realizations of the same
process and exchanging the subsystems. However, there are cases where it is not
possible to measure several realizations, like, e.g., in geophysical systems.
The construction of the surrogates we present in this section is also based on
the recurrence matrix (10.4). It is important to note that if the recurrence matrix
is computed from a univariate time series, it contains all topological information
about the underlying attractor, which therefore can be reconstructed from it [29].
Hence, a first idea for the generation of surrogates is to change the structures
in a RP consistently with the ones produced by the underlying dynamical system.
In this way one could reconstruct a new realization of the trajectory from the
modified Ri,j. However, one cannot arbitrarily interchange columns in an RP,
because such a modification changes the distribution of diagonal lines and hence
the entropy and predictability of the system [28].
Therefore, we propose a modified approach. In general, in an RP there are
identical columns, i.e., Rk,i = Rk,j ∀k [28]. Thus, there are points which are not
only neighbors (i.e., ‖xi − xj‖∞ < ε), but which also share the same neighbor-
hood. Reconstructing the attractor from an RP, the respective neighborhoods of
these points cannot help to distinguish them, i.e., from this point of view they
are identical. This is why we will call them twins. Twins are special points of the
time series as they are indistinguishable considering their neighborhoods but in
general different and hence, have different pasts and—more important—different
futures. The key idea of how to introduce the randomness needed for the gener-
ation of surrogates of a deterministic system is that one can jump randomly to
one of the possible futures of the existing twins.

266 10 Synchronization Analysis and Recurrence in Complex Systems
A surrogate trajectory xs(i) of x(i) with i = 1, . . . , N is then generated in the
following way:
1. Identify all pairs of twins.
2. Choose an arbitrary starting point, say xs(1) = x(k).
3. If x(k) has no twin, the next point of the surrogate trajectory is xs(2) = x(k+1).
4. If x(k) has a twin, say x(m), then one can go to either x(k + 1) or to x(m + 1),
i.e., xs(2) = x(k + 1) or xs(2) = x(m + 1) with equal probability6.
Steps three and four are then iterated until the surrogate time series has the same
length as the original one.
This algorithm creates twin surrogates (TS) which are shadows of a (typical)
trajectory of the system [64]. In the limit of an infinitely long original trajectory,
its surrogates are characterized by the same dynamical invariants and the same
attractor. However, if the measure of the attractor can be estimated from the ob-
served finite trajectory reasonably well, its surrogates share the same statistics.
Also their power spectra and correlation functions are consistent with the ones
of the original system. TS do not only seem to give reasonable results for deter-
ministic systems; the TS of for example an ARMA process also show the typical
behavior of a linear Gaussian process.
Next, we use the TS to test for PS. The idea behind this approach is simi-
lar to the one by means of “natural surrogates” in the mother–fetus heartbeat
synchronization [15]. Suppose that we have two coupled self-sustained oscilla-
tors x1(t) and x2(t). Then, we generate M TS of the joint system, i.e., xsi1 (t)
and xsi2 (t), with i = 1, . . . , M. These surrogates are independent copies of the
joint system, i.e., trajectories of the whole system beginning at different initial
conditions. Note that the coupling between x1(t) and x2(t) is also mimicked
by the surrogates. Next, we compute the differences between the phases of the
original system ∆Φ(t) = Φ1(t) − Φ2(t) applying, e.g., the analytical signal ap-
proach [1–4] and compare them with ∆Φsi(t) = Φ1(t)−Φ
si
2 (t) (one can also con-
sider Φsi1 (t) − Φ2(t)). Then, if ∆Φ(t) does not differ significantly from ∆Φ
si(t)
with respect to some index for PS, the null hypothesis cannot be rejected and
hence, we do not have enough evidence to state PS.
As a test case, we consider two nonidentical, mutually coupled Rössler oscil-
lators
x˙1,2 = −(1 ± ν)y1,2 − z1,2 + ε(x2,1 − x1,2),
y˙1,2 = (1 ± ν)x1 + 0.15y1,2,
z˙1,2 = 0.2 + z1,2 + z1,2(x1,2 − 10) ,
(10.14)
where ν = 0.015 denotes the frequency mismatch. In this “active experiment”, we
vary the coupling strength ε from 0 to 0.08 and compute a PS index for the orig-
inal trajectory for each value of ε. Next we generate 200 TS and compute the PS
6 If triplets occur one proceeds analogously.

10.5 Twin Surrogates to Test for PS 267
index between the measured first oscillator and the surrogates of the second one.
As PS index we use the mean resultant length R of complex phase vectors [65, 66],
which is motivated by Kuramoto’s order parameter [67]
R =
∣
∣
∣
∣
∣
1
N
N
∑
t=1
exp
(
i∆Φ(t)
)
∣
∣
∣
∣
∣
. (10.15)
It takes on values in the interval from 0 (non PS) to 1 (perfect PS) [65, 66]. Let Rsi
denote the PS index between the first oscillator and the surrogate i of the second
one. To reject the null hypothesis at a significance value α, R must be larger
than (1 − α) · 100 percent of all Rsi . Note that this corresponds to computing the
significance level from the cumulative histogram at the level (1 − α).
Figure 10.18(a) shows the results for R of the original system (bold line) and
the 1% (solid) significance level. Figure 10.18(b) displays the difference between R
of the original system and the 1%, 2% and 5% significance level. For ε < 0.025, R
of the original system is, as expected, below the significance levels and hence the
difference is negative, and for higher values of ε the curves cross (the difference
becomes positive). This is in agreement with the criterion for PS via Lyapunov
exponents λi [1–4]: λ4 becomes negative at ε ∼ 0.028 (Fig. 10.18(b)), which ap-
proximately coincides with the intersection of the curve of R for the original
system and the significance level (zero-crossing of the curves in Fig. 10.18(b)).
Therefore, we recognize successfully the PS region by means of the TS.
Note that also the significance limit increases when the transition to PS occurs
(Fig. 10.18(a)). As the TS mimic both the linear and nonlinear characteristics of
the system, the surrogates of the second oscillator have in the PS region the same
mean frequency as the first original oscillator. Hence Rsi is rather high. However,
Φ1(t) and Φ
si
2 (t) do not adapt to each other, as they are independent. Hence, the
value of R for the original system is significantly higher than the Rsi . We state
in conclusion that even though the obtained value for a normalized PS index
is higher than 0.97 (right side of Fig. 10.18(a)), this does not offer conclusive
evidence for PS. Hence, the knowledge of the PS index alone does not provide sufficient
evidence for PS. Note that the more phase coherent the oscillators are, the more
difficult it is to decide whether they are in PS or not. A certain phase diffusion,
which allows to measure the adaptation of the phases of the interacting oscillators
is necessary to detect PS. However, the test based on the TS reveals whether there
is enough evidence for PS.
Next, we perform an analysis of the specificity and sensitivity of the test
for ε = 0 and ν = 0. For 100 random initial conditions of the Rössler system
and a significance level of α = 1%, the null hypothesis was erroneously rejected
only in 1 out of the 100 cases. This is a rather auspicious result, as due to the
identical frequencies, it is extremely difficult to recognize that there is no PS in
this case [68]. In the case of ε = 0.02 (e.g., no PS) and ν = 0.015, there were
no erroneous rejections of the null hypothesis. Finally, for PS (ε = 0.045 and

268 10 Synchronization Analysis and Recurrence in Complex Systems
Fig. 10.18: (a) R of the original two mutually coupled Rössler systems with a fre-
quency mismatch of ν = 0.015 (bold) and significance level of 1% (solid). (b) Dif-
ference between R of the original data and significance level of 1% (solid), 2%
(dashed) and 5% (dashed-dotted). The zero line is plotted (dotted) to guide the
eye. (c) Four largest Lyapunov exponents for the six-dimensional system consid-
ered. λ4 is highlighted and the arrow indicates the transition to PS.
ν = 0.015), in all 100 test runs the null hypothesis was correctly rejected. These
results indicate that the specificity and sensitivity of the test are good.
10.6 Application to Fixational Eye Movements
Next we apply the recurrence approach to check fixational movements of left and
right eyes for PS. During fixation of a stationary target our eyes perform small
involuntary and allegedly erratic movements to counteract retinal adaptation. If
these eye movements are experimentally suppressed, retinal adaptation to the
constant input induces very rapid perceptual fading [69, 70]. Moreover, statistical
correlations show a timescale separation from persistence to antipersistence [71].
Persistence on the short timescale counteracts retinal fading, whereas antipersis-
tence on the long timescale contributes to stability of ocular disparity. According
to current textbook knowledge, the fixational movements of the left and right eyes
are correlated very poorly at best [72]. Therefore, it is highly desirable to exam-
ine these processes from a perspective of PS. We analyze the data of two subjects.
Each performed three trials, in which they fixated a small stimulus (black square
on a white background, 3× pixels on a computer display) with a spatial extent of

10.6 Application to Fixational Eye Movements 269
Fig. 10.19: Simultaneous recording of left (bold) and right (solid) fixational eye
movements (a) horizontal component (b) vertical component.
0.12°, or 7.2 arc ·min. Eye movements were recorded using an EyeLink-II system
(SR Research, Osgoode, Ontario, Canada) with a sampling rate of 500Hz and
an instrument spatial resolution less than 0.005°. Figure 10.19 shows a segment
of the horizontal (a) and vertical (b) component of the eye movements for one
person.
The data were first high-pass filtered applying a difference filter x˜(t) = x(t) −
x(t − τ) with τ = 40ms in order to eliminate the slow drift of the data. Af-
ter this filtering, we find an oscillatory trajectory, which has maximum spectral
power in the frequency range between 3 and 8Hz (Fig. 10.20(a) and (b)). How-
ever, the trajectories of the eyes are rather noisy and non-phase-coherent. There-
fore, it is cumbersome to estimate the phase of these data. Hence we apply the
recurrence-based measure CPR introduced in Section 10.2 and we obtain the val-
ues displayed in the first column of Table 10.2. First, we observe that the variabil-
ity between the different trials is smaller for the first participant as for the second
one. Furthermore, the values of CPR are rather high for the first participant but
not so high for the second one. Hence, a hypothesis test should be performed in
order to get statistically significant results.
Therefore, we compute 200 twin surrogates of the left eye’s trajectory. In
Fig. 10.20(c) the horizontal component of one surrogate is represented. At a first
glance, the characteristics of the original time series are well reproduced by the
twin surrogate. In Fig. 10.20(d) the corresponding periodogram is displayed. It
is also noteworthy that the structure of the original curve (Fig. 10.20(b)) is also
qualitatively reproduced. The periodogram of the twin surrogate is of course not
identical with the one of the original time series. This is consistent with the null
hypothesis of another realization of the same underlying process, respectively
another trajectory starting at different initial conditions of the same underlying
dynamical system.
Now, we compute the recurrence-based synchronization index CPRsi between
the twin surrogates of the left eye and the measured right eye’s trajectory. In
Fig. 10.21 the results of the test of one trial are visualized.

270 10 Synchronization Analysis and Recurrence in Complex Systems
Fig. 10.20: Filtered horizontal component of the left eye of one participant (a) and
its corresponding periodogram (b). In (c) the horizontal component of one surro-
gate of the left eye is represented and in (d) its corresponding periodogram.
Fig. 10.21: Result of the test performed for one trial of one participant. The PS
index for the original data (bold line) is significantly different from the one of the
surrogates (solid).
The second column of Table 10.2 summarizes the results for both subjects and
all trials.
In all cases, the PS index of the original data is significantly different from the
ones of the surrogates, which strongly indicates that the concept of PS can be suc-
cessfully applied to study the interaction between the trajectories of the left and

10.7 Conclusions 271
Tab. 10.2: Results for the test for PS between the trajectories of the left and right
fixational eye movements performed for three trials for the two participants. Two
hundred surrogates were used for the test. The null hypothesis was rejected in all
cases at a 2% level.
Participant CPR of the original data Null hypothesis
M.R. 0.9112 Rejected
0.9432 Rejected
0.9264 Rejected
M.T. 0.6080 Rejected
0.4844 Rejected
0.3520 Rejected
right eyes during fixation. This result also suggests that the physiological mech-
anism in the brain that produces the fixational eye movements controls both eyes
simultaneously, i.e., there might be only one center in the brain that produces
the fixational movements in both eyes or a close link between two centers. Our
finding of PS between left and right eyes is in good agreement with current
knowledge of the physiology of the oculomotor circuitry. In a single-cell study,
66% of abducens motor neurons fired in relation to the movements of either eye,
while premotor neurons in the brainstem encode monocular movements [73].
Thus, motor neurons – as the final common pathway of neural control of eye
movements – are candidates for the synchronization of left and right fixational
movements. Furthermore, we are interested in whether the fixational movements
in the horizontal and vertical direction of one eye are synchronized. Horizon-
tal and vertical saccadic eye movements are controlled in two spatially distinct
brainstem nuclei [74]. Therefore, we can expect that, on the level of fixational
eye movements, horizontal and vertical components are independent. Applying
the synchronization index CPR between the x- and y-component of the left eye of
each participant for each trial and generate 200 surrogates of the two-dimensional
trajectory of the left eye. Then we compare the synchronization index CPRsi be-
tween the original x-component and the y-component of the surrogates. We find
in all but one cases that CPR is not significantly different from CPRsi (see Ta-
ble 10.3). Hence, we do not have evidence to claim synchronization between the
horizontal and vertical components of the eye movements, as expected.
10.7 Conclusions
In conclusion, we have presented solutions to two main problems of the syn-
chronization analysis of measured time series: The detection of PS in non-phase-
coherent systems and the hypothesis test for PS, which is interesting especially
for passive experiments, where the coupling strength between the two subsys-
tems cannot be varied systematically.

272 10 Synchronization Analysis and Recurrence in Complex Systems
Tab. 10.3: Results for the test for PS between the horizontal and vertical compo-
nents of fixational movements of one eye performed for three trials for the two
participants. 200 surrogates were used for the test. In all cases but one, we failed
to reject the null hypothesis at a 2% level.
Participant CPR of the original data Null hypothesis
M.R. 0.3746 Not rejected
0.6103 Not rejected
0.4812 Rejected
M.T. 0.4681 Not rejected
0.3194 Not rejected
0.4172 Not rejected
We have given solutions to these two problems based on the concept of recur-
rence in phase space. First, we have shown that by means of the recurrence prop-
erties it is possible to detect indirectly PS even in the case of non-phase-coherent
and strong noisy time series. Furthermore, it is also possible to detect GS by
means of recurrences. Second, the method of twin surrogates has been presented,
which is also based on recurrence, and we have shown that it can be used to test
for PS.
We have used the well studied system of two mutually coupled Rössler oscil-
lators in order to validate the techniques proposed. Furthermore, we have tested
for PS in experiments of binocular fixational movements and found that the left
and right eyes are in PS, in agreement with physiological results about the func-
tional role of motor neurons in the final common pathway for the control of eye
movements. Hence, we have shown that the techniques proposed are also ap-
plicable to rather noisy observed time series.
Acknowledgements
This work has been suported by the DFG Priority Program 1114 and the “Interna-
tionales Promotionskolleg—Helmholtz Center for the Study of Mind and Brain
Dynamics” at the University of Potsdam.Author: Please pro-
vide publisher and
place in [3]
Author: Please up-
date [63]
References
[1] M. Rosenblum, A. Pikovsky, and J. Kurths. Phys. Rev. Lett., 76:1804, 1996.
[2] A. Pikovsky, M. Rosenblum, G. Osipov, and J. Kurths. Physica D, 104:219,
1997.
[3] A. Pikovsky, M. Rosenblum, and J. Kurths. Synchronization, volume 12 of
Cambridge Nonlinear Science Series. 2001.

10.7 Conclusions 273
[4] S. Boccaletti, J. Kurths, G. V. Osipov, D. Valladares, and C. Zhou. Phys. Rep.,
366:1, 2002.
[5] R. C. Elson et al. Phys. Rev. Lett., 81:5692, 1998.
[6] P. Tass et al. Phys. Rev. Lett., 81:3291, 1998.
[7] C. M. Ticos et al. Phys. Rev. Lett., 85:2929, 2000.
[8] V. Makarenko and R. Llinas. Proc. Natl. Acad. Sci. USA, 95:15474, 1998.
[9] B. Blasius, A. Huppert, and L. Stone. Nature, 399:354, 1999.
[10] C. Schäfer et al. Nature, 392:239, 1998.
[11] D. J. DeShazer et al. Phys. Rev. Lett., 87:044101, 2001.
[12] S. Boccaletti et al. Phys. Rev. Lett., 89:194101, 2002.
[13] I. Kiss and J. Hudson. Phys. Rev. E, 64:046215, 2001.
[14] H. Poincaré. Acta Math., 13:1, 1890.
[15] P. Van Leeuwen et al. BMC Physiol., 3:2, 2003.
[16] M. G. Rosenblum, A. S. Pikovsky, J. Kurths, G. V. Osipov, I. Z. Kiss, and J. L.
Hudson. Phys. Rev. Lett., 89:264102, 2002.
[17] G. V. Osipov, B. Hu, C. Zhou, M. V. Ivanchenko, and J. Kurths. Phys. Rev.
Lett., 91:024101, 2003.
[18] G. Fisher. Plane Algebraic Curves. American Mathematical Society, Provi-
dence, Rhode Island, 2001.
[19] C. Sparrow. The Lorenz Equations: Bifurcations, Chaos, and Strange Attractors.
Springer, Berlin, Heidelberg, 1982.
[20] R. N. Madan. Chua Circuit: A Paradigm for Chaos. World Scientific, Singapore,
1993.
[21] W. Lauterborn, T. Kurz, and M. Wiesenfeldt. Coherent Optics. Fundamentals
and Applications. Springer, Berlin, Heidelberg, New York, 1993.
[22] J. Y. Chen et al. Phys. Rev. E, 64:016212, 2001.
[23] J.-P. Eckmann, S. O. Kamphorst, and D. Ruelle. Europhys. Lett., 4:973, 1987.
[24] C. L. Weber Jr. and J. P. Zbilut. J. Appl. Physiol., 76:965, 1994.
[25] N. Marwan, N. Wessel, U. Meyerfeldt, A. Schirdewan, and J. Kurths. Phys.
Rev. E, 66:026702, 2002.

274 10 Synchronization Analysis and Recurrence in Complex Systems
[26] N. Marwan and J. Kurths. Phys. Lett. A, 302:299, 2002.
[27] M. Thiel et al. Physica D, 171:138, 2002.
[28] M. Thiel, M. C. Romano, P. Read, and J. Kurths. Chaos, 14:234, 2004.
[29] M. Thiel, M. C. Romano, and J. Kurths. Phys. Lett. A, 330:343, 2004.
[30] E.-H. Park, M. Zaks, and J. Kurths. Phys. Rev. E, 60:6627, 1999.
[31] V. S. Afraimovich, N. N. Verichev, and M. I. Rabinovich. Izvestiya Vysshikh
Uchebnykh Zavedenii Radiofizika, 29:1050, 1986.
[32] N. F. Rulkov, M. M. Sushchik, L. S. Tsimring, and H. D. I. Abarbanel. Phys.
Rev. E, 51 2:980, 1995.
[33] U. Parlitz and L. Kocarev. Synchronization of chaotic systems. In H. G.
Schuster, editor, Handbook of Chaos Control. Wiley-VCH, Weinheim, 1999.
[34] L. M. Pecora and T. L. Carroll. Phys. Rev. Lett., 64:821, 1990.
[35] K. M. Cuomo and A. V. Oppenheim. Phy. Rev. Lett., 71:65, 1993.
[36] A. Kittel, A. Parisi, and K. Pyragas. Physica D, 112:459, 1998.
[37] G. D. Van Wiggeren and R. Roy. Science, 279:1198, 1998.
[38] K. Otsuka, R. Kawai, S.-L. Hwong, J.-Y. Ko, and J.-L. Chern. Phys. Rev. Lett.,
84:3049, 2000.
[39] C. W. Wu and L. O. Chua. Int. J. Bif. Chaos, 4:1979, 1994.
[40] T. L. Carroll and L. M. Pecora. Physica D, 67:126, 1993.
[41] L. Kocarev and U. Parlitz. Phys. Rev. Lett., 74:5028, 1995.
[42] U. Parlitz, L. Kocarev, T. Stojanovski, and H. Preckel. Phys. Rev. E, 53:4351,
1996.
[43] L. Kocarev, K. S. Halle, K. Eckert, L. O. Chua, and U. Parlitz. Int. J. Bif. Chaos,
2:709, 1992.
[44] T. L. Carroll and L. M. Pecora. Physica D, 67:126, 1993.
[45] U. Parlitz, L. Junge, and L. Kocarev. Phys. Rev. E, 54:6253, 1996.
[46] J. Arnhold, P. Grassberger, K. Lehnertz, and C. E. Elger. Physica D, 134:419,
1999.
[47] R. Q. Quiroga, J. Arnhold, and P. Grassberger. Phys. Rev. E, 61:5142, 2000.
[48] A. Schmitz. Phys. Rev. E, 62:7508, 2000.

10.7 Conclusions 275
[49] S. J. Schiff, P. So, T. Chang, R. E. Burke, and T. Sauer. Phys. Rev. E, 54:6708,
1996.
[50] M. Wiesenfeldt, U. Parlitz, and W. Lauterborn. Int. J. Bif. Chaos, 11:2217,
2001.
[51] M. C. Romano, M. Thiel, J. Kurths, and W. von Bloh. Phys. Lett. A, 330:214,
2004.
[52] M. G. Rosenblum, A. S. Pikovsky, and J. Kurths. Phys. Rev. Lett., 78:4193,
1997.
[53] L. Kocarev and U. Parlitz. Phys. Rev. Lett., 76:1816, 1996.
[54] F. Takens. Detecting strange attractors in turbulence. In D. A. Rand and
L.-S. Young, editors, Dynamical Systems and Turbulence, volume 898 of Lecture
Notes in Mathematics. Springer, Berlin, 1980.
[55] H. Kantz and T. Schreiber. Nonlinear Time Series Analysis. Cambridge Uni-
versity Press, Cambridge, 1997.
[56] M. C. Romano, M. Thiel, J. Kurths, I. Z. Kiss, and J. Hudson. Europhys. Lett.,
71:466, 2005.
[57] D. Maraun and J. Kurths. Geophys. Res. Lett., 32:15709, 2005.
[58] M. Palus. Phys. Lett. A, 235:341, 1997.
[59] M. Palus and A. Stefanovska. Phys. Rev. E, 67:055201(R), 2003.
[60] J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, and J. D. Farmer. Physica D,
58:77, 1992.
[61] T. Schreiber and A. Schmitz. Phys. Rev. Lett., 77:635, 1996.
[62] M. Small, D. Yu, and R. G. Harrison. Phys. Rev. Lett., 87:188101, 2001.
[63] M. Thiel, M. C. Romano, J. Kurths, and R. Engbert. submitted.
[64] E. Ott. Chaos in Dynamical Systems. Cambridge University Press, Cambridge,
1993.
[65] E. Rodriguez et al. Nature, 397:430, 1999.
[66] C. Allefeld and J. Kurths. Int. J. Bif. Chaos, 14:405, 2004.
[67] Y. Kuramoto. Chemical Oscillations, Waves and Turbulence. Springer, New
York, 1984.
[68] M. Peifer, B. Schelter, M. Winterhalder, and J. Timmer. Phys. Rev. E, 72:
026213, 2005.

276 10 Synchronization Analysis and Recurrence in Complex Systems
[69] L. A. Riggs, F. Ratliff, J. C. Cornsweet, and T. N. Cornsweet. J. Opt. Soc. Am.,
43:495, 1953.
[70] D. Coppola and D. Purves. Proc. Natl. Acad. Sci. USA, 93:8001, 1996.
[71] R. Engbert and R. Kliegl. Psychol. Science, 15:431, 2004.
[72] K. J. Ciuffreda and B. Tannen. Eye Movement Basics for the Clinician. Mosby,
St. Louis, 1995.
[73] W. Zhou and W. M. King. Nature, 393:692, 1998.
[74] D. L. Sparks. Nature Rev. Neurosci., 3:952, 2002.

Index
L1-regression, 154, 156
α-trimmed mean filter, 148
AIC, 392
Akaike Information Criterion (AIC),
462
Akaike information criterion (AIC), 392
Akaike’s Final Prediction Error (FPE),
392
analytic signal, 24, 26, 176, 228
AR models, 91
ARMA(X) Systems, 301
artificial neural networks, 84
attractor, 82
Autoregressive model, 282
autoregressive model, 452, 461
Autoregressive representation, 349
autoregressive-moving average model,
391
Bayesian Information Criterion (BIC),
463
Bayesian theorem, 96
Bias Variance Dilemma, 53
biosurveillance, 322
biosurveillance data, 332
bivariate data, 171, 178
block bootstrap, 434
Brain Machine Interface, 128
breakdown point, 144, 146–148, 150,
154–156, 160
Burg (LWR) algorithm, 393
cardiorespiratory coordination, 174
cardiorespiratory interaction, 178
causal influence, 453, 456
causality, 185
CCA-Subspace Estimators, 307
chaotic oscillators, 247, 264
Circuit Data, 285
climatic data, 213
Cluster weighted modeling, 65
coherence, 393, 455
conditional Granger causality, 459, 461
connectivity, 427
Contemporaneous correlation, 350
continuity measure, 281
correlation abolishing TR transforma-
tions, 430
correlation of probability of recurrence,
253
correlation preserving TP transforma-
tions, 430
coupled chaotic oscillators, 232
coupled oscillators, 171, 172, 175, 212
Cross Validation, 55
cross validation, 41, 43, 45
cross-correlation, 278
cross-correlation analysis, 172
cross-correlation function, 172, 189
cross-spectral analysis, 172
curvature, 247
data driven modeling, 295
dDTF, 395
delay embedding, 23, 30, 44
delay in coupling, 172, 187
depth electroencephalography, 234
Detection of coupling, 277
Determinism, 101
determinism, 81
DFT – Discrete Fourier Transform, 432
diagnostics of coupling, 212
direct Directed Transfer Function
(dDTF), 395
Directed coherence (DC), 427
Directed transfer function, 371
– direct — (dDTF), 373
Directed transfer function (DTF), 427
directed transfer function (DTF), 395
direction of coupling, 212
directionality index, 186, 187

506 Index
directionality of coupling, 172
directionality of interaction, 183, 185
double window, 165
double window filter, 149, 150, 152, 153,
156, 157
DTF, 395
Dynamic Linear Model (DLM), 324
dynamical systems, 81, 193
edge, 147–149, 152, 158, 159, 161–163
EEG, 475, 478
effective brain connectivity, 475
El Niño/Southern Oscillation, 214
electrocardiogram, 179
electrocardiograms, 101
electroencephalogram, 215
electroencephalograms, 101
EM Algorithm, 67
Embedding, 112
embedding, 82
EMD applications, 226
entropy measures, 185
epilepsy, 234
epilepy, 215
Task Force of the European Society of
Cardiology, 140
events, 93
exact fit point, 146, 147, 154
exponentially weighted moving aver-
age, 145
factor models, 296
Factor Models for Time Series, 311
Factor Models with Idiosyncratic Noise,
313
false nearest neighbors, 17, 24, 25, 27
fixational eye movements, 268
fMRI, 475
Fokker-Planck equation, 88
FPE, 392
frequency mismatch, 185
Generalized Linear Dynamic Factor
Models, 315
Generalized PDC (GPDC), 427
generalized synchronization, 256
global embedding dimension, 24
Global Principal Component Analysis,
110, 114
Granger Causality, 475
Granger causality, 349–354, 394, 395,
426, 451
– bivariate —, 353
– multivariate —, 350
Granger causality graph, 355
– bivariate —, 357
Granger’s causality concept, 185
Hénon map, 284
high dimensional time series, 295
high-dimensional time series, 101
Hilbert transform, 24, 26, 176, 228, 281
Hilbert-Huang transform, 227
horizon of predictability, 17
hybrid filter, 143, 157, 159, 160, 162, 165
IDFT – Inverse Discrete Fourier Trans-
form, 433
iid – independent and identically
distributed, 428
impulse detection, 161, 163
instantaneous causality, 454, 457
instantaneous phase, 228
intensity of interaction, 183
interaction, 172
interdependence, 453, 455
interspike intervals, 30, 40
intrinsic time scales, 227
invertible model, 431
joint probability of recurrence, 258
Kalman filter, 325, 396
Langevin equation, 88
Latent variable, 365
least median of squares, 154, 164
level 2 statistics, 428
level shift, 144–147, 149, 150, 154–158,
160, 161, 163, 164
limit cycle, 174
linear filter, 143, 145, 147, 152, 390
linear least-square regression, 180
local embedding dimension, 17
Local Modeling, 51
LWR, 393

Index 507
Lyapunov exponent, 175
Lyapunov exponents, 18, 31, 35, 36
m-separation, 359
magnetoencephalography, 235
Markov chain, 89
Markov process, 89
Markov property, 363, 364
mathematical modeling, 193, 218
MDL, 117
melanoma incidence, 440
membrane voltage, 18, 19, 21, 23, 41, 44
mesial temporal lobe epilepsy, 442
Minimum Description Length, 117
modeling, 17, 40, 41
models, global nonlinear, 84
modified trimmed means filter, 149, 155
modulation, 175, 178
Morlet, or Gabor, wavelet, 177
moving average, 143, 145, 149, 150, 390
moving window, 145, 156
multichannel measurements, 87
multiple coherence, 394
multiple shooting, 45
multivariate autoregressive model, 425
multivariate autoregressive models, 390
mutual entrainment, 174
Mutual Information, 279, 280
mutual information, 17
mutual predictability, 185
nearest neighbors prediction, 17
neural synchronization, 226
neuron time series, 40
Noise, 104
non-phase-coherent oscillators, 248
nonautonomous systems, 217
normalized Directed Transfer Function,
395
North Atlantic Oscillation, 214
online, 145–147, 153, 155, 156, 165
order statistic filter, 148, 149, 151, 165
outlier, 143–147, 149, 154–156, 158, 160,
164
Overfitting, 53
overfitting, 40, 42, 43, 45
oversampling, 81
parametric models, 387
partial coherence, 393
Partial coherence (PC), 427
Partial directed coherence, 353, 367
Partial directed coherence (PDC), 427
partial directed coherence (PDC), 395
Partial directed correlation, 367
Partial spectral coherency, 360
Path diagram, 354
– bivariate —, 357
PDC, 395
periodically forced systems, 231
permutation procedure, 472
phase, 174
phase and frequency locking, 174, 178
Phase correlation, 280
phase diffusion, 264
phase dynamics, 175
phase resampling, 432
phase shift, 172
phase slips, 247
phase synchronization, 175, 261, 264,
268, 280
phase-locking index, 233
Poincaré section, 184
point process, 177
predictability improvement, 186
Prediction, 52
prediction error, 186
prediction errors, 85, 91
prediction: Markov chain, 90
predictions, more step, 86
predictor, locally constant, 83
predictor, locally linear, 83
Principal Component Analysis, 312
Principal Component Regression, 60
probabilistic prediction, 95
probability of recurrence, 249
Rössler in funnel regime, 252
Rössler system, 284
radial basis functions, 42
Randomness, 102, 104
randomness, 101
reconstruction, 193, 206
recurrence plot, 248
recursive filter, 145, 148, 159

508 Index
regularization, 43
reliability test, 95
repeated median, 154–157, 159, 163–165
repeated median filter, 155–157, 159,
165
residue resampling, 431
respiratory sinus arrhythmia, 178
Ridge regression, 60
robust filter, 143, 144, 156
robust regression, 143, 153, 155, 156,
163, 165
ROC statistics, 94, 96
root signal, 148, 152, 159, 164
running median, 143, 145–152, 158, 159,
161, 162, 164
sampling rate, 81
second order statistics, 434
seizure focus, 442
self-sustained oscillator, 172, 173
semi-nonparametric identification, 305
shift detection, 161–163
shift-dependent synchronization index,
188
signal extraction, 143, 146, 153, 164
signal processing, 101
Spatial Granger Causality, 480
Spatially Constrained Models, 328
spectral distribution function, 297
spectral matrix, 390, 463
Spectral representation, 352
speech, 101
spike, 143–147, 149, 151, 153, 158,
160–163
STARMAX model, 321, 329, 331
State Space Model, 324
state space reconstruction, 30
state space systems, 301
stationary processes, 296
Stochastic, 277
Stochastic resonance, 106
strange attractor, 174
strength of coupling, 184
stroboscopic approach, 184
stroboscopic synchronization index, 184
surrogate hypothesis testing, 182
synchrogram, 184
synchronization, 174
synchronization index, 183, 188
Synchrony, 277
system identification, 193
time delay embedding, 82
time scale synchronization, 232
time series, 193
transfer matrix, 390
trend, 143–148, 150, 153, 155, 157–161,
163–165
twin surrogates, 265
Uncertainty, 105
update algorithm, 155, 156, 159, 165
Vector autoregressive model
– graphical —, 369
Visual Evoked Potentials, 131
Wölfer sunspot data, 440
weather prediction, 79
weighted median filter, 151, 152, 156
Wold decomposition, 298
Yule–Walker algorithm, 392