diff --git a/.gitignore b/.gitignore
index 175ec13..eb0e6a2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,4 +4,4 @@ Manifest.toml
 *.css
 *style.jl
 docs/src/examples/*.md
-paper/wip
\ No newline at end of file
+paper/data/*.csv
\ No newline at end of file
diff --git a/ext/TransitionVisualizations.jl b/ext/TransitionVisualizations.jl
index dad3ad7..f263260 100644
--- a/ext/TransitionVisualizations.jl
+++ b/ext/TransitionVisualizations.jl
@@ -51,8 +51,6 @@ function TransitionsInTimeseries.plot_indicator_changes(res::SegmentedWindowResu
     fig, n = init_rowwise_visualisation(res, colors, indicator_names,
         chametric_names, accent_linewidth)
     for k in eachindex(res.t_indicator)
-        Makie.lines!(contents(fig[1, 1])[1], res.t[k], res.x[k], color = colors[1],
-            linewidth = accent_linewidth)
         lineplot_metrics!(fig, n, res.config, res.t_indicator[k], res.x_indicator[k],
             res.t_change[k], res.x_change[k, :], colors, accent_linewidth)
     end
@@ -65,17 +63,17 @@ function init_rowwise_visualisation(res, colors, indicator_names,
     chametric_names, accent_linewidth)
 
     fig = Makie.Figure(size = (700, 450), fontsize = 12)
-    rlabels = vcat([""], indicator_names)
-    llabels = vcat(["Input"], chametric_names)
+    llabels = vcat(["Input"], indicator_names)
+    rlabels = vcat([""], chametric_names)
     n = length(llabels)
 
     # rowaspect = 5
-    raxs = [Makie.Axis(fig[i, 1], ylabel = rlabels[i],
-        xticklabelsvisible = false, yaxisposition = :right, ygridvisible = false,
-        ylabelcolor = colors[2], yticklabelcolor = colors[2]) for i in 1:n]
     laxs = [Makie.Axis(fig[i, 1], ylabel = llabels[i],
         xticklabelsvisible = false, ygridvisible = false,
         ylabelcolor = colors[1], yticklabelcolor = colors[1]) for i in 1:n]
+    raxs = [Makie.Axis(fig[i, 1], ylabel = rlabels[i],
+        xticklabelsvisible = false, yaxisposition = :right, ygridvisible = false,
+        ylabelcolor = colors[2], yticklabelcolor = colors[2]) for i in 1:n]
     Makie.linkxaxes!(laxs..., raxs...)
 
     hidedecorations!(raxs[1])
@@ -90,7 +88,7 @@ function init_rowwise_visualisation(res, colors, indicator_names,
     labels = ["original signal", "surro signals"]
     width = 0.5
     if length(res.t_indicator[1]) > 1
-        elements = vcat(elements, [MarkerElement(marker = :circle, color = colors[1],
+        elements = vcat(elements, [MarkerElement(marker = :circle, color = colors[2],
             strokecolor = :transparent, markersize = ms) for ms in [10, 5]])
         labels = vcat(labels, ["original change metric", "surro change metrics"])
         width = 1
@@ -106,14 +104,14 @@ function lineplot_metrics!(fig, n, config, t_ind, x_ind, t_cha, x_cha,
     for i in 2:n
         j = i-1
         if !isnothing(config.indicators)
-            Makie.lines!(contents(fig[i, 1])[2], t_ind, x_ind[:, j], color = colors[2],
+            Makie.lines!(contents(fig[i, 1])[1], t_ind, x_ind[:, j], color = colors[1],
                 linewidth = accent_linewidth)
         end
         if length(t_cha) > 1
-            lines!(contents(fig[i, 1])[1], t_cha, x_cha[:, j], color = colors[1],
+            lines!(contents(fig[i, 1])[2], t_cha, x_cha[:, j], color = colors[2],
                 linewidth = accent_linewidth)
         else
-            Makie.scatter!(contents(fig[i, 1])[1], t_cha, x_cha[j], color = colors[1],
+            Makie.scatter!(contents(fig[i, 1])[2], t_cha, x_cha[j], color = colors[2],
                 markersize = 10)
         end
     end
@@ -166,26 +164,28 @@ function lines_over_surro!(fig, nsurro, t, t_ind, t_cha, x, signif, config,
     for ns in 1:nsurro
         s = TimeseriesSurrogates.surrogate(x, signif.surrogate)
         Makie.lines!(contents(fig[1, 1])[1], t, s; color = (colors[1], 2/nsurro), linewidth = 1)
-        for (i, cha) in enumerate(config.change_metrics)
 
+        for (i, cha) in enumerate(config.change_metrics)
             if isnothing(config.indicators)
                 p = s
             else
                 p = windowmap(config.indicators[i], s; width = config.width_ind,
                     stride = config.stride_ind)
-                Makie.lines!(contents(fig[i+1, 1])[2], t_ind, p; color = (colors[2], 2/nsurro),
+                Makie.lines!(contents(fig[i+1, 1])[1], t_ind, p; color = (colors[1], 2/nsurro),
                     linewidth = 1)
             end
+
             if length(t_cha) > 1
                 q = windowmap(cha, p; width = config.width_cha, stride = config.stride_cha)
-                Makie.lines!(contents(fig[i+1, 1])[1], t_cha, q; color = (colors[1], 2/nsurro),
+                Makie.lines!(contents(fig[i+1, 1])[2], t_cha, q; color = (colors[2], 2/nsurro),
                     linewidth = 1)
             else
                 cha = precompute(cha, eachindex(p))
                 q = windowmap(cha, p; width = length(p), stride = 1)
-                Makie.scatter!(contents(fig[i+1, 1])[1], t_cha, q[1], color = (colors[1], 2/nsurro),
+                Makie.scatter!(contents(fig[i+1, 1])[2], t_cha, q[1], color = (colors[2], 2/nsurro),
                     markersize = 5)
             end
+            
             if !isnothing(flags) && ns == 1 && length(t_cha) > 1
                 Makie.vlines!(contents(fig[i+1, 1])[1], t_cha[flags[:, i]];
                     color = (:black, 0.1))
diff --git a/paper/code/Project.toml b/paper/code/Project.toml
index 9e5445f..59cf384 100644
--- a/paper/code/Project.toml
+++ b/paper/code/Project.toml
@@ -1,4 +1,4 @@
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
-TransitionsInTimeseries = "5f5b98ec-1183-43e0-887a-12fdc55c52f7"
diff --git a/paper/code/ewstools-tuto-1.py b/paper/code/ewstools-tuto-1.py
index 457fba8..f65c0b5 100644
--- a/paper/code/ewstools-tuto-1.py
+++ b/paper/code/ewstools-tuto-1.py
@@ -3,58 +3,74 @@
 import matplotlib.pyplot as plt
 import ewstools
 from ewstools.models import simulate_ricker
+from importlib.metadata import version
 
-# Set seed for reproducibility
-np.random.seed(0)
+if version('ewstools') != '2.1.0':
+    raise ValueError('Please install version 2.1.0 of ewstools')
 
-# Initialize time series and spectrum computation
-series = simulate_ricker(tmax=1000, F=[0,2.7])
-ts = ewstools.TimeSeries(data=series, transition=860)
-ts.detrend(method='Lowess', span=0.2)
-ts.state[['state','smoothing']].plot()
-ts.compute_spectrum(rolling_window=0.5, ham_length=40)
+def ewstools_setup():
+    # Initialize time series and spectrum computation
+    series = simulate_ricker(tmax=1000, F=[0,2.7])
+    ts = ewstools.TimeSeries(data=series, transition=860)
+    ts.detrend(method='Lowess', span=0.2)
+    ts.state[['state','smoothing']].plot()
+    ts.compute_spectrum(rolling_window=0.5, ham_length=40)
+    return ts
 
-# Initialize parameters for timing functions
-rw = 0.5
-n = 100
-t_elapsed = np.zeros(10)
+np.random.seed(0)           # Set random seed for reproducibility
+rw = 0.5                    # Rolling window for variance computation
+ts = ewstools_setup()       # Setup time series
+n = 100                     # Number of iterations for timing
+t_elapsed = np.zeros(n)     # Array to store elapsed times
+t_minruntime = np.zeros(8) # Array to store minimum runtime
 
-# Time functions (in a not very elegant way)
-t0 = time.time()
 for i in range(n):
+    t0 = time.time()
     ts.compute_var(rolling_window=rw)
-t_elapsed[0] = time.time() - t0
+    t_elapsed[i] = time.time() - t0
+t_minruntime[0] = min(t_elapsed)
 
-t0 = time.time()
 for i in range(n):
+    t0 = time.time()
     ts.compute_cv()
-t_elapsed[1] = time.time() - t0
+    t_elapsed[i] = time.time() - t0
+t_minruntime[1] = min(t_elapsed)
 
-t0 = time.time()
 for i in range(n):
+    t0 = time.time()
     ts.compute_skew(rolling_window=rw)
-t_elapsed[2] = time.time() - t0
+    t_elapsed[i] = time.time() - t0
+t_minruntime[2] = min(t_elapsed)
 
-t0 = time.time()
 for i in range(n):
+    t0 = time.time()
     ts.compute_kurt()
-t_elapsed[3] = time.time() - t0
+    t_elapsed[i] = time.time() - t0
+t_minruntime[3] = min(t_elapsed)
 
-t0 = time.time()
 for i in range(n):
+    t0 = time.time()
     ts.compute_auto(rolling_window=rw, lag=1)
-t_elapsed[4] = time.time() - t0
+    t_elapsed[i] = time.time() - t0
+t_minruntime[4] = min(t_elapsed)
 
-t0 = time.time()
 for i in range(n):
+    t0 = time.time()
     ts.compute_smax()
-t_elapsed[5] = time.time() - t0
+    t_elapsed[i] = time.time() - t0
+t_minruntime[5] = min(t_elapsed)
 
-t0 = time.time()
 for i in range(n):
+    t0 = time.time()
     ts.compute_ktau()
-t_elapsed[6] = time.time() - t0
+    t_elapsed[i] = time.time() - t0
+t_minruntime[6] = min(t_elapsed)
 
-t0 = time.time()
-surro = ewstools.core.block_bootstrap(ts.state.residuals, n, bs_type='Stationary', block_size=10)
-t_elapsed[7] = time.time() - t0
\ No newline at end of file
+for i in range(n):
+    t0 = time.time()
+    ewstools.core.block_bootstrap(ts.state.residuals, 1, bs_type='Stationary', block_size=10)
+    t_elapsed[i] = time.time() - t0
+t_minruntime[7] = min(t_elapsed)
+
+np.savetxt('../data/ewstools_ricker.csv', ts.state[['residuals']].values)
+np.savetxt('../data/ewstools_perfo.csv', t_minruntime, delimiter=',')
\ No newline at end of file
diff --git a/paper/code/figure1.jl b/paper/code/figure1.jl
index 148fd3e..0168dc2 100644
--- a/paper/code/figure1.jl
+++ b/paper/code/figure1.jl
@@ -1,6 +1,8 @@
 using TransitionsInTimeseries, DelimitedFiles, CairoMakie, Random
 
-x = readdlm("ewstools-tuto-1.csv", ',')[:, end]
+# Run ewstools-tuto-1.csv first to generate ewstools_ricker.csv
+x = readdlm("../data/ewstools_ricker.csv", ',')[:, end]
+x = x[isnan.(x) .== 0]
 t = eachindex(x)
 
 # Choose the indicators and how to measure their change over time
@@ -12,5 +14,4 @@ results = estimate_changes(config, x, t)
 signif = SurrogatesSignificance(n = 1000, tail = [:right, :right], rng = Xoshiro(1995))
 flags = significant_transitions(results, signif)
 fig = plot_changes_significance(results, signif)
-ylims!(contents(fig[2, 1])[2], (0.037, 0.045))
 save("../figures/figure1.png", fig)
\ No newline at end of file
diff --git a/paper/code/figure2.jl b/paper/code/figure2.jl
index 74be5f5..3e7e427 100644
--- a/paper/code/figure2.jl
+++ b/paper/code/figure2.jl
@@ -1,69 +1,60 @@
-using TransitionsInTimeseries, DelimitedFiles, CairoMakie, Random
+using CairoMakie
+using DelimitedFiles
+using Random
+using TransitionsInTimeseries
 
 coefficient_of_variation(x) = std(x) / mean(x)
 
 function main()
-    x = readdlm("ewstools-tuto-1.csv", ',')[:, end]
+    # Run ewstools-tuto-1.csv first to generate ewstools_ricker.csv and ewstools_perfo.csv
+    x = readdlm("../data/ewstools_ricker.csv", ',')[:, end]
+    x = x[isnan.(x) .== 0]
     t = eachindex(x)
 
     # Choose the indicators and how to measure their change over time
     indicators = (var, coefficient_of_variation, skewness, kurtosis,
         ar1_whitenoise, LowfreqPowerSpectrum())
     stride = [1, 1, 1, 1, 1, 40]
-    n, m = 100, length(indicators)
+    m = length(indicators)
     t_elapsed = zeros(m+2)
 
-
     for (i, ind) in enumerate(indicators)
         # Build configuration with adequate parameters of the sliding window
         config = SegmentedWindowConfig((ind, ind), (nothing, nothing), [t[1]], [t[end]];
             width_ind = length(x) ÷ 2, stride_ind = stride[i], whichtime = last,
             min_width_cha = 1)
 
-        t0 = time()
-        for i in 1:n
-            # Compute the metrics over sliding windows and their significance
-            results = estimate_changes(config, x, t)
-        end
-        t_elapsed[i] = (time() - t0) / 2
+        t_elapsed[i] = @belapsed estimate_changes($config, $x, $t)
     end
 
     config = SegmentedWindowConfig((nothing, nothing), (kendalltau, kendalltau), [t[1]], [t[end]];
         width_ind = length(x) ÷ 2, stride_ind = 1, whichtime = last, min_width_cha = 1)
-    t0 = time()
-    for i in 1:n
-        results = estimate_changes(config, x, t)
-    end
-    t_elapsed[m+1] = (time() - t0) / 2
+    t_elapsed[m+1] = @belapsed estimate_changes($config, $x, $t)
 
     sgen = surrogenerator(x, BlockShuffle(), Xoshiro(1995))
-    t0 = time()
-    for i in 1:n
-        s = sgen()
-    end
-    t_elapsed[m+2] = time() - t0
+    t_elapsed[m+2] = @belapsed $sgen()
 
     return t_elapsed
 end
 
-t_tt = main()
-t_et = [0.03840542, 0.05554581, 0.03895116, 0.04029274, 7.96556187,
-    2.73067856, 0.39529872, 0.02751493]
-
-# [0.04681492, 8.13679838, 0.04035759, 0.09219241]
-inds = eachindex(t_et)
+t_transitionsintimeries = main()
+t_ewstools = readdlm("../data/ewstools_perfo.csv", ',')[:, end]
+inds = eachindex(t_ewstools)
 w = 0.4
 
-fig, ax = barplot(inds .- 0.5*w, t_et, label = L"ewstools $\,$", width = w,
-    fillto = 1e-5)
-barplot!(ax, inds .+ 0.5*w, t_tt, label = L"TransitionsInTimeseries.jl $\,$",
-    width = w, fillto = 1e-5)
+fig = Figure(resolution = (600, 400))
+ax = Axis(fig[1, 1])
+ylims!(ax, 1e-7, 0.1)
+barplot!(ax, inds .- 0.5*w, t_ewstools, label = L"ewstools $\,$", width = w,
+    fillto = 1e-8)
+barplot!(ax, inds .+ 0.5*w, t_transitionsintimeries, label = L"TransitionsInTimeseries.jl $\,$",
+    width = w, fillto = 1e-8)
 ax.yscale = log10
 ax.xticks = (1:8, [L"Variance $\,$", L"Coeff. of variation $\,$", L"Skewness $\,$",
     L"Kurtosis $\,$", L"Lag-1 autocorr. $\,$", L"Spectral $\,$",
     L"Kendall $\tau$ corr. coeff.", L"Block bootstrap $\,$"])
-ax.ylabel = L"Run time (s) of 100 computations on Ricker model data $\,$"
-ax.yticks = (10.0 .^ (-5:1), [L"10^{%$e}" for e in -5:1])
+ax.ylabel = L"Run time (s) on Ricker model data $\,$"
+ax.yticks = (10.0 .^ (-8:1), [L"10^{%$e}" for e in -8:1])
 ax.xgridvisible = false
 ax.ygridvisible = false
 ax.xticklabelrotation = π / 4
diff --git a/paper/figures/figure1.png b/paper/figures/figure1.png
index cf84ae8..9f9bfd8 100644
Binary files a/paper/figures/figure1.png and b/paper/figures/figure1.png differ
diff --git a/paper/figures/figure2.png b/paper/figures/figure2.png
index 93def1c..8443246 100644
Binary files a/paper/figures/figure2.png and b/paper/figures/figure2.png differ
diff --git a/paper/paper.md b/paper/paper.md
index 4ea106d..ee099ac 100644
--- a/paper/paper.md
+++ b/paper/paper.md
@@ -44,12 +44,15 @@ work predicting a collapse of the Atlantic Meridional Overturning Circulation be
 and 2095 has led to no less than 870 news outlets and 4100 tweets [@ditlevsen_warning_2023],
 largely because of the substantial implications of such a collapse for human societies.
 A common concern in the scientific community is that published work on the topic is difficult
-to reproduce, despite the impact it implies for humanity.
-This can be largely addressed by a unifying software that is accessible, performant,
-reproducible, reliable and extensible. Such a software does not exist yet, but here
-we propose TransitionsInTimeseries.jl to fill this gap.
-We believe this is a major step towards establishing a software as standard, widely used
-by academics working on transitions in timeseries.
+to reproduce, despite the impact it implies for humanity. Additionally, existing softwares
+[@dakos_methods_2012, @bury_ewstools_2023] are well suited to apply some state-of-the-art
+techniques but lack the possibility of being easily extended by the users, thus making them
+unsuited for some research task.
+Both of these issues can be largely addressed by a unifying software that is accessible,
+performant, reproducible, reliable and extensible. Such a software does not exist yet, but
+here we propose TransitionsInTimeseries.jl to fill this gap. We believe this is a major step
+towards establishing a software as standard, widely used by academics working on transitions
+in timeseries.
 
 # TransitionsInTimeseries.jl
 
@@ -77,8 +80,11 @@ surrogates, corresponding to a p-value $p<0.05$. All these steps can be performe
 visualisation of the results within a few lines only:
 
 ```julia
+using TransitionsInTimeseries, CairoMakie, Random
+
 # Loading and preprocessing the data needs to be done by the user
 time, data = load_data()
+residual = data - filt(data)
 
 # Choose the indicators and how to measure their change over time
 indicators = (var, ar1_whitenoise)
@@ -89,8 +95,8 @@ config = SegmentedWindowConfig(indicators, change_metrics, [time[1]], [time[end]
     width_ind = length(residual) ÷ 2, whichtime = last, min_width_cha = 100)
 
 # Compute the metrics over sliding windows and their significance
-results = estimate_changes(config, data, time)
-signif = SurrogatesSignificance(n = 1000, tail = :right, rng = Xoshiro(1995))
+results = estimate_changes(config, residual, time)
+signif = SurrogatesSignificance(n = 1000, tail = (:right, :right), rng = Xoshiro(1995))
 flags = significant_transitions(results, signif)
 
 # Visualize the results
@@ -118,7 +124,7 @@ TransitionsInTimeseries.jl is written in Julia, which offers both a simple synta
 performance. Additionally, all performance-relevant steps have been optimized and
 parallelized when possible, as, for instance, the significance testing relying on surrogates.
 In the final section of this article, we present a comparison to `ewstools`, showing that
-TransitionsInTimeseries.jl offers a significant speed-up in all the studied cases.
+TransitionsInTimeseries.jl offers a significant speed-up in most the studied cases.
 
 ## Reproducibility
 
@@ -173,17 +179,15 @@ thus offering optimized routines with numerous surrogate types.
 ### Choosing a pipeline
 
 TransitionsInTimeseries.jl covers methods for prediction as well as detection of transitions,
-which is unprecedented to our knowledge. This relies on the definition of different analysis pipelines, which
-consist in a `ChangesConfig` determining the behavior of `estimate_changes` via
-multiple dispatch. For instance, a detection
-task can be performed by replacing the `SegmentedWindowConfig` by a `SlidingWindowConfig`
-in the code above:
+which is unprecedented to our knowledge. This relies on the definition of different analysis
+pipelines, which consist in a `ChangesConfig` determining the behavior of `estimate_changes`
+via multiple dispatch. For instance, a detection task can be performed by replacing the `SegmentedWindowConfig` by a `SlidingWindowConfig` in the code above:
 
 ```julia
 # Here the data should not be detrended
 time, data = load_data()
 
-indicators = (nothing, nothing)
+indicators = nothing
 change_metrics = (difference_of_mean(), difference_of_max())
 config = SlidingWindowConfig(indicators, change_metrics;
    width_cha = 50, whichtime = midpoint)
@@ -210,10 +214,9 @@ significance pipeline and makes TransitionsInTimeseries.jl particularly versatil
 # Comparison to already existing alternatives
 
 `earlywarnings` [@dakos_methods_2012] and `spatialwarnings` are toolboxes written in R
-providing many tools to predict transitions.
-These are early and valuable efforts but are (1) restricted to
-prediction tasks, (2) written in a less performant language, (3) not parallelised, (4) not
-designed for convenient reproducibility and (5) not extensible.
+providing many tools to predict transitions. These are early and valuable efforts but
+are (1) restricted to prediction tasks, (2) written in a less performant language,
+(3) not parallelised, (4) not designed for convenient reproducibility and (5) not extensible.
 
 `ewstools` [@bury_ewstools_2023] is a Python/TensorFlow package offering similar
 functionalities as `earlywarnings`, as well as a deep-learning approach to predicting
@@ -224,11 +227,11 @@ detection and prediction tasks. We believe that this is now covered by
 TransitionsInTimeseries.jl.
 
 Using TransitionsInTimeseries.jl, we reproduced the computations showcased in Tutorial 1
-and Tutorial 2 of the `ewstools` documentation, along with the block bootstrapping. We
-performed each computation 100 times and show the resulting run times in [Fig. 2](@figure2).
-It appears that all computation are faster in TransitionsInTimeseries, with a speed-up factor
-ranging from 0 to 3 orders of magnitude. The implementation of the deep-learning classifiers
-for transition prediction developed in [@bury_deep_2021], as well as dealing with
+and Tutorial 2 of the `ewstools` (v2.1.0) documentation, along with the block bootstrapping.
+The runtimes of both softwares are benchmarked in [Fig. 2](@figure2).
+It appears that most computations are faster in TransitionsInTimeseries, with a speed-up
+factor ranging from 0 to 3.5 orders of magnitude. The implementation of the deep-learning
+classifiers for transition prediction developed in [@bury_deep_2021], as well as dealing with
 multidimensional timeseries, are part of future developments of TransitionsInTimeseries.jl.
 
 ![Performance comparison between `ewstools` and TransitionsInTimeseries.jl.\label{fig: fig1}](figures/figure2.png)
@@ -241,8 +244,8 @@ The documentation of TransitionsInTimseries.jl is available at
 # Acknowledgements
 
 Jan Swierczek-Jereczek is funded by CriticalEarth, grant no. 956170, an H2020 Research
-Infrastructure of the European Commission. George Datseris is funded by UKRI's Engineering and
-Physical Sciences Research Council, grant no. EP/Y01653X/1 (grant agreement for a EU
+Infrastructure of the European Commission. George Datseris is funded by UKRI's Engineering
+and Physical Sciences Research Council, grant no. EP/Y01653X/1 (grant agreement for a EU
 Marie Sklodowska-Curie Postdoctoral Fellowship).
 
 # References
\ No newline at end of file
diff --git a/src/misc/timeseries.jl b/src/misc/timeseries.jl
index 2e0b023..d248f32 100644
--- a/src/misc/timeseries.jl
+++ b/src/misc/timeseries.jl
@@ -7,10 +7,10 @@ differences in `t` are compared with `isapprox(; kwargs...)` and if all are appr
 the same then `true` is returned.
 """
 function isequispaced(t::AbstractVector; kwargs...)
-    # TODO: assert 1 based indexing
-    s = t[2] - t[1]
-    for i in 3:length(t)
-        sn = t[i] - t[i-1]
+    fi = firstindex(t)
+    s = t[fi+1] - t[fi]
+    for i in 2:length(t)-1
+        sn = t[fi+i] - t[fi+i-1]
         isapprox(sn, s; kwargs...) || return false
     end
     return true
@@ -19,9 +19,10 @@ isequispaced(t::AbstractRange; reltol = 1e-6) = true
 
 equispaced_step(t::AbstractRange) = step(t)
 function equispaced_step(t::AbstractVector)
+    fi = firstindex(t)
     s = zero(eltype(t))
-    for i in 2:length(t)
-        s += t[i] - t[i-1]
+    for i in 1:length(t)-1
+        s += t[fi+i] - t[fi+i-1]
     end
     return s/(length(t)-1)
 end
diff --git a/test/Project.toml b/test/Project.toml
index ec1a036..870bb6e 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,11 +1,9 @@
 [deps]
-CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
-DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
-Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TimeseriesSurrogates = "c804724b-8c18-5caa-8579-6025a0767c70"
-StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
\ No newline at end of file
diff --git a/test/timeseries.jl b/test/timeseries.jl
index 942edc8..2554e5c 100644
--- a/test/timeseries.jl
+++ b/test/timeseries.jl
@@ -1,12 +1,14 @@
-using TransitionsInTimeseries, Test
+using TransitionsInTimeseries, Test, OffsetArrays
 
 @testset "spacing" begin
     dt = rand()
     x_equispaced = range( 0.0, step = dt, stop = 100)
+    x_equispaced_offset = OffsetVector([1:10...], -1:8)
     x_nonequispaced = cumsum( rand(1000) )
     @test isequispaced(x_equispaced)
     @test !isequispaced(x_nonequispaced)
     @test equispaced_step(x_equispaced) ≈ dt
+    @test equispaced_step(x_equispaced_offset) ≈ 1
 end
 
 #=