S_map_Sugihara1994 {pttstability}

Apply S-mapping algorithm from Sugihara 1994

Description

Carries out an S-mapping analysis, following the algorithm outlined in Sugihara (1994).

Usage

S_map_Sugihara1994(Y, E, theta, X = NULL, lib = NULL, trimNA = FALSE)

Arguments

Y

a timeseries vector from which to build the embedding.

E

a positive integer, specifying the embedding dimension

theta

a positive numeric scalar, specifying the nonlinearity parameter for the analysis. A value of 0 indicates a fully linear analysis; higher numbers indicate greater nonlinearity.

X

an optional matrix of time-delayed embeddings to use for the analysis

lib

an optional matrix of library positions, for specifying cases where Y is a composite timeseries made up of multiple separate observations (e.g. spatial replicates). Matrix should have two columns, with the first row in each column specifying the start of the timeseries section, and the second column specifying the end.

trimNA

a logical specifying whether NA values should be removed from Y and X - defaults to FALSE

Value

a list, including the timeseries used for S-mapping (Y), the delay embedding matrix used for S-mapping (X), a vector of predictions (Y_hat), a matrix of S-mapping coefficients (C), the standard errors for the S-mapping coefficients (C_SE), and goodness of fit metrics R-squared (R2) and root mean square error (RMSE).

Source

Sugihara, G. (1994). Nonlinear forecasting for the classification of natural time-series. Philos. Trans. R. Soc. -Math. Phys. Eng. Sci., 348, 477–495.

Examples

# create an example timeseries
n = 100
set.seed(1234)
datout<-makedynamics_general(n = n+2,
                             pdet=log(c(0.8,1)),
                             proc = -2.5,
                             detfun = detfun0_sin)
plot(datout$true, type = "l") # plot timeseries
Y = datout$true # extract true values

# run s-mapping
sout = S_map_Sugihara1994(Y = Y, E = 2, theta = 0.5)
s_coef = process_scof(sout$C) # process coefficients from the S-mapping output

# find best E/theta
fitout = data.frame(E = 1:5, theta = NA, RMSE = NA)

for(i in 1:nrow(fitout)) {
   E = fitout$E[i]
   Ytmp = Y[-c(1:E)]
   optout = optimize(f = function(x) {S_map_Sugihara1994(Ytmp, E, x)$RMSE}, interval = c(0,10))
  
  fitout$theta[i] = optout$minimum # get best theta for given E
  fitout$RMSE[i] = optout$objective # get error
}
ps = which.min(fitout$RMSE)

E = fitout$E[ps] # get best E
theta = fitout$theta[ps] # get best theta
X = makeblock(Y, E) # get X for analysis
Y = Y[-c(1:E)] # trim NA values (corresponding to positions in X)
X = X[(E+1):nrow(X),] # trim NA values
sout = S_map_Sugihara1994(Y = Y, E = E,
  theta = theta, X = X) # run S-mapping for best paramter combination
sout$R2 # look at R-squared

# check fit
plot(sout$Y_hat, Y)
abline(a=0, b=1, lty=2)