Questions

• Use lx <- lorenz[1:2048,1] to reconstruct the phase space based on lx.
  • Find an optimal embedding lag and dimension by calling crqa_parameters().
  • Use nnThres = .01. This is the proportion of remaining nearest neighbours below which we consider the number of dimensions optimal. Because this time series was generated by a simulation of a system, we can use a much lower threshold. For real data, which is much noisier, you can use 0.1 or 0.2.
  • NOTE: If you have a slow computer don’t use 2048, but use e.g. 512.

This crqa_parameters() function will call several other functions (e.g. fractal::timeLag() with method = "mutual" and nonlinearTseries::findAllNeighbours() for the false nearest neighbourhood search). Look at the manual page, you can run the crqa_parameters() function with the defaults settings, it will show a diagnostic plot. Because we have to search over a range of parameters it can take a while to get the results…

• Based on the output:
  • Choose an embedding lag \(\tau\) (in casnet functions this is the argument emLag). Remember, this is just an optimisation choice. Sometimes there is no first minimum, or no global minimum.
  • Find an appropriate embedding dimension by looking when the false nearest neighbours drop below 1%.
• Now you can use this information to reconstruct the phase space:
  • Embed the time series using ts_embed() (you can also use fractal::embedSeries()).
  • Inspect the object that was returned
  • Use rgl::plot3d() to plot the reconstructed space. Plot the reconstructed phase space. (You’ll need to use as.matrix() on the object created by fractal::embedSeries())

Answers

• Use lx <- lorenz[1:2048,1] to reconstruct the phase space based on lx.
  • Find an optimal embedding lag and dimension by calling crqa_parameters().

library(rgl)
library(ggplot2)
library(fractal)

## Loading required package: splus2R

## Loading required package: ifultools

library(invctr)
library(casnet)

# The X data
lx <- lorenz[1:2048,1]

# Search for the parameters,
params <- crqa_parameters(lx)

## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo

# Assign the optimal parameters
emLag <- params$optimLag
emDim <- params$optimDim

• Now you can use this information to reconstruct the phase space

# For X
lx_emb <- ts_embed(y=lx,emLag = emLag, emDim = emDim)

rgl::plot3d(lx_emb, type="l")

# For Y
ly <- lorenz[1:2048,2]

# Search for the parameters,
params <- crqa_parameters(ly)

# Assign the optimal parameters
emLag <- params$optimLag
emDim <- params$optimDim

# Embed the time series
ly_emb <- ts_embed(y=ly,emLag = emLag, emDim = emDim)
rgl::plot3d(ly_emb,type="l")

# For Z
lz <- lorenz[1:2048,3]

# Search for the parameters,
params <- crqa_parameters(lz)

# Assign the optimal parameters
emLag <- params$optimLag
emDim <- params$optimDim

# Embed the time series
lz_emb <- ts_embed(y=lz,emLag = emLag, emDim = emDim)
rgl::plot3d(lz_emb,type="l")

Questions

• First get a radius that gives you a fixed recurrence rate use crqa_radius() look at the manual pages.

• Best way to ensure you are using the same parameters in each function is to create some lists with parameter settings (check the crqa_cl() manual to figure out what these parameters mean)

Answers

library(fractal)
library(casnet)

# Lorenz X
lx <- lorenz[1:2048,1]
emLag = 17
emDim = 3

(emRad <- crqa_radius(y1 = lx, emLag = emLag, emDim = emDim)$Radius)

## 
## Auto-recurrence: Setting diagonal to (1 + max. distance) for analyses

## lower and upper are both 0 (no band, just diagonal)
##  using: diag(mat) <- 38.5536...

## 
## Searching for a radius that will yield 0.05 for RR

## 
## Converged! Found an appropriate radius...

## [1] 3.629007

# RQA analysis
(out <- crqa_cl(lx,emDim = emDim, emLag = emLag, emRad = emRad))

## 
## ~~~o~~o~~casnet~~o~~o~~~
## 
## Performing auto-RQA
## 
## ...using sequential processing...
## 
  |                                                                     
  |                                                               |   0%
  |                                                                     
  |o~oo~oo~oo~oo~oo~oo~oo~oo~oo~oo~oo~oo~oo~oo~oo~oo~oo~oo~oo~oo~o| 100%
## 
## Completed in:
##    user  system elapsed 
##   1.972   0.020   1.998 
## 
## ~~~o~~o~~casnet~~o~~o~~~

##                                 .id        RR      DET  DET.RR     LAM
## 1 window: 1 | start: 1 | stop: 2048 0.0499637 0.998963 19.9938 0.99918
##   LAM.DET L_max       L  L_entr         DIV V_max      TT  V_entr      T1
## 1 1.00022  2013 24.6951 3.76081 0.000496771    35 10.1375 2.75883 15.9085
##       T2 W_max  W_mean  W_entr W_prob       F_min emDim emLag    emRad
## 1 175.25  1533 164.849 4.67858     58 0.000652316     3    17 3.629007
##   DLmin VLmin  distNorm
## 1     2     2 EUCLIDEAN

We can plot the recurrence matrix…

library(casnet)

# Unthresholded matrix (no radius applied)
RM <- rp(y1 = lx, emDim = emDim, emLag = emLag)

# plot it
rp_plot(RM)

# Thresholded by the radius
RMth <- di2bi(RM,emRad = emRad)

rp_plot(RMth, plotDimensions = TRUE, plotMeasures = TRUE)

Questions

Use RQA to compare white noise to the deterministic chaos series of previous assignments.

library(rio)
library(pracma)

# Reload the data
series <- rio::import("https://github.com/FredHasselman/The-Complex-Systems-Approach-Book/raw/master/assignments/assignment_data/BasicTSA_arma/series.sav", setclass = "tbl_df")

Answers

Discussed in the next session

Questions

• Perform an RQA on one of the coordinate series (x or y)
• Study what happens to the RQA measures if you shuffle the temporal order.
• Package fractal contains a function surrogate. This will create so-called constrained realisations of the time series. You can look at the vignette in package casnet cl_CRQA (also available in Github) Look at the help pages of the function, or study the Surrogates Manual of the TISEAN software and create two surrogate series, one based on phase and one on aaft.
  • Look at the RQA measures and think about which \(H_0\) should probably be rejected.
  • If you want to be more certain, you’ll have to create a test (more surrogates). The TISEAN manual provides all the info you need to construct such a test:

  “For a minimal significance requirement of 95% , we thus need at least 19 or 39 surrogate time series for one- and two-sided tests, respectively. The conditions for rank based tests with more samples can be easily worked out. Using more surrogates can increase the discrimination power.”