Some months ago, I posted a function that extended the base R function `embed()` to allow for time delay embedding. Today, David Gonzales alerted me to an inconsistency between `embed()` and `Embed()`.

The example David used was

```R> (x <- seq(1,20,3))
  1  4  7 10 13 16 19
R> embed(x, 4)
[,1] [,2] [,3] [,4]
[1,]   10    7    4    1
[2,]   13   10    7    4
[3,]   16   13   10    7
[4,]   19   16   13   10
R> Embed(x, 4)
[,1] [,2] [,3] [,4]
[1,]    4    3    2    1
[2,]    7    6    5    4
[3,]   10    9    8    7
[4,]   13   12   11   10
```

where `Embed()` clearly returns an incorrect result.

In this post, I present an explanation of the problem and address the shortcomings in the original code with an updated version of `Embed()`.

The reason the original version of `Embed()` doesn’t work with David’s example is that when I wrote it, I had in mind that it would work on the indices of the time series, not the values of the time series. I had overlooked that `embed()` returned the embedded time series, not the indices — the problem of testing with vectors like `1:10`!

Updating `Embed()` to output the same result as `embed()` is a trivial matter; we just get the function to work with `seq_along(x)` and not `x` itself and then use the old `Embed()` behaviour to index `x` to return the embedded time series. As an added extra, as we are generating the indices anyway, we can optionally have the function return those instead of the embedded series.

Here is the updated version of `Embed()`

```Embed <- function(x, m, d = 1, indices = FALSE, as.embed = TRUE) {
n <- length(x) - (m-1)*d
X <- seq_along(x)
if(n <= 0)
stop("Insufficient observations for the requested embedding")
out <- matrix(rep(X[seq_len(n)], m), ncol = m)
out[,-1] <- out[,-1, drop = FALSE] +
rep(seq_len(m - 1) * d, each = nrow(out))
if(as.embed)
out <- out[, rev(seq_len(ncol(out)))]
if(!indices)
out <- matrix(x[out], ncol = m)
out
}
```

The main difference is that we create `X <- seq_along(x)` and create `out` using that rather than the time series (`x`). I’ve also added a new argument, `indices`, that defaults to `FALSE`. If we want `Embed()` to return the indices of the embedded time series, call the function with `indices = FALSE`.

The new version of `Embed()` gives the same results as before and is consistent with `embed()` when we pass it a time series that is identical to its indices

```R> embed(1:5, 2)
[,1] [,2]
[1,]    2    1
[2,]    3    2
[3,]    4    3
[4,]    5    4
R> Embed(1:5, 2)
[,1] [,2]
[1,]    2    1
[2,]    3    2
[3,]    4    3
[4,]    5    4
```

but it also works for time series like those in David’s example:

```R> (x <- seq(1,20,3))
  1  4  7 10 13 16 19
R> embed(x, 4)
[,1] [,2] [,3] [,4]
[1,]   10    7    4    1
[2,]   13   10    7    4
[3,]   16   13   10    7
[4,]   19   16   13   10
R> Embed(x, 4)
[,1] [,2] [,3] [,4]
[1,]   10    7    4    1
[2,]   13   10    7    4
[3,]   16   13   10    7
[4,]   19   16   13   10
```

and we have the added benefit of being able to return the indices of the embedded time series

```R> Embed(x, 4, indices = TRUE)
[,1] [,2] [,3] [,4]
[1,]    4    3    2    1
[2,]    5    4    3    2
[3,]    6    5    4    3
[4,]    7    6    5    4
```

Now I just need to do something on the recurrence plot that I originally wrote `Embed()` for!