# Numerically estimate the Jacobian from a scalar time-series

I'm trying to numerically estimate the Jacobian from a time-series. Following the paper, https://journals.aps.org/pra/abstract/10.1103/PhysRevA.43.2787.

Consider that I have a scalar time series $$x = (x(1), x(2), \dots, x(N - 1), x(N))$$. In total there are $$N$$ data-points. So the steps that I'm following are as follows.

1. Compute the time-delay (call it $$\tau$$).
2. Compute the embedding dimension (call it $$d_{emb}$$).
3. Compute the delay vectors (call it $$\bar y_{i}(k)$$). $$\bar y_{i}(k) = [x(i), x(i + \tau), x(i + 2\tau), x(i + 3\tau), \dots, x(i + (d_{emb} - 1) \tau]$$.
4. Choose some random initial point in the time series and call the delay vector corresponding to that point, $$\bar y(k)$$.
5. Look at the neighbours of $$\bar y(k)$$ and store them in ascending order of their Euclidean distance in $$\bar y^{(r)}(k)$$. Say we keep $$N_{b}$$ neighbours.
6. Store the distance between the delay vector and its neighbours in $$\bar z^{(r)}(k) = \bar y^{(r)}(k) - \bar y(k)$$.

The basic idea is that the underlying dynamics is governed by a map $$\bar F$$ which takes $$\bar y(k) \rightarrow \bar y(k + 1)$$.

1. Compute the next evolution of $$\bar y(k)$$. Therefore, $$\bar y(k+1) = \bar F(\bar y(k))$$. Similarly do it for the $$r$$th neighbour of $$\bar y(k)$$. One gets, $$\bar y(r; k+1) = \bar F(\bar y^{(r)}(k))$$.
2. Compute the distance between these vectors, $$\bar z(r; k+1) = \bar y(r; k+1) - \bar y(k+1)$$.
3. Now, $$\bar z(r; k+1) = \bar y(r; k+1) - \bar y(k + 1) = \bar F(y^{(r)}(k)) - \bar F(y(k)) = \bar F(\bar y(k) + \bar z^{(r)}(k)) - \bar F(\bar y(k))= M_{\alpha \beta} z^{(r)}_{\beta}(k) + Q_{\alpha \beta \gamma} z^{(r)}_{\beta}(k) z^{(r)}_{\gamma}(k) + \cdots$$.

We assume that $$\bar z^{(r)}(k)$$ are small so that one can do a Taylor expansion about it. $$M$$ is the Jacobian matrix of dimension, $$((d_{emb} - 1)\tau, (d_{emb} - 1)\tau)$$. $$Q$$ is the Hessian matrix.

1. Idea now is to compute the unknown coefficients of $$M$$ using least-square fit keeping $$H$$ (for accuracy and dropping it finally).

The authors suggest that if $$\tau$$ be the step-size for the next iteration, so $$\bar y(k) \rightarrow \bar y(k + 1)$$, then $$\bar z(r; k + 1) = M_{\alpha \beta} z^{(r)}(k)$$. And then the only components unknown are the bottom row of $$M$$ Rest are $$1$$'s along the superdiag, and rest are all $$0$$.

What I don't understand it is that $$\bar z^{(r)}(k) \rightarrow \bar F(\bar z^{(r)}(k)) = \bar F(\bar y^{(r)}(k) - \bar y(k)) \ne \bar F(\bar y^{(r)}(k)) - \bar F(\bar y (k)) = z(r; k+1)$$. So how do they write -- $$\bar z(r; k+1) = M_{\alpha \beta} z^{(r)}_{\beta}(k)$$?