diff --git a/.gitignore b/.gitignore
index b22a5e4d3d..1143cec0ee 100644
--- a/.gitignore
+++ b/.gitignore
@@ -165,3 +165,7 @@ local_code/
 local/
 local_code.py
 local.py
+
+/benchmarks/benchmark_results/
+
+benchmarks/html/
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 847bcf11b7..d68ef8bd26 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -31,6 +31,7 @@ repos:
     hooks:
       - id: ruff
         args: [ "--fix"]
+        exclude: '(^|/)benchmarks/'
 
   - repo: https://github.com/asottile/pyupgrade
     rev: v3.20.0
diff --git a/benchmarks/README.md b/benchmarks/README.md
new file mode 100644
index 0000000000..40607a227b
--- /dev/null
+++ b/benchmarks/README.md
@@ -0,0 +1,74 @@
+[//]: # (This was adapted from: https://github.com/scipy/scipy/tree/main/benchmarks)
+# aeon Time Series Benchmarks
+
+Benchmarking aeon with Airspeed Velocity.
+
+## Usage
+
+Airspeed Velocity manages building and Python environments by itself, unless told
+otherwise. Some of the benchmarking features in `spin` also tell ASV to use the aeon
+compiled by `spin`. To run the benchmarks, you will need to install the "dev"
+dependencies of aeon:
+
+```bash
+pip install --editble .[dev]
+# NOTE: If the above fails, try running pip install --editable ".[dev]"
+```
+
+Run a benchmark against currently checked-out aeon version (don't record the result):
+
+```bash
+spin bench --submodule classification.distance_based
+```
+
+Compare change in benchmark results with another branch:
+
+```bash
+spin bench --compare main --submodule classification.distance_based
+```
+
+Run ASV commands directly (note, this will not set env vars for `ccache` and disabling BLAS/LAPACK multi-threading, as `spin` does):
+
+```bash
+cd benchmarks
+asv run --skip-existing-commits --steps 10 ALL
+asv publish
+asv preview
+```
+
+More on how to use `asv` can be found in [ASV documentation](https://asv.readthedocs.io/). Command-line help is available as usual via `asv --help` and `asv run --help`.
+
+## Writing benchmarks
+
+See [ASV documentation](https://asv.readthedocs.io/) for the basics on how to write benchmarks.
+
+Some things to consider:
+
+- When importing things from aeon on the top of the test files, do it as:
+
+  ```python
+  from .common import safe_import
+
+  with safe_import():
+      from aeon.classification.distance_based import KNeighborsTimeSeriesClassifier
+  ```
+
+  The benchmark files need to be importable also when benchmarking old versions of aeon. The benchmarks themselves don't need any guarding against missing features — only the top-level imports.
+
+- Try to keep the runtime of the benchmark reasonable.
+
+- Use ASV's `time_` methods for benchmarking times rather than cooking up time measurements via `time.clock`, even if it requires some juggling when writing the benchmark.
+
+- Preparing arrays etc., should generally be put in the `setup` method rather than the `time_` methods, to avoid counting preparation time together with the time of the benchmarked operation.
+
+- Use `run_monitored` from `common.py` if you need to measure memory usage.
+
+- Benchmark versioning: by default `asv` invalidates old results when there is any code change in the benchmark routine or in setup/setup_cache.
+
+  This can be controlled manually by setting a fixed benchmark version number, using the `version` attribute. See [ASV documentation](https://asv.readthedocs.io/) for details.
+
+  If set manually, the value needs to be changed manually when old results should be invalidated. In case you want to preserve previous benchmark results when the benchmark did not previously have a manual `version` attribute, the automatically computed default values can be found in `results/benchmark.json`.
+
+- Benchmark attributes such as `params` and `param_names` must be the same regardless of whether some features are available, or e.g. AEON_XSLOW=1 is set.
+
+  Instead, benchmarks that should not be run can be skipped by raising `NotImplementedError` in `setup()`.
diff --git a/benchmarks/asv.conf.json b/benchmarks/asv.conf.json
new file mode 100644
index 0000000000..a02fad59af
--- /dev/null
+++ b/benchmarks/asv.conf.json
@@ -0,0 +1,87 @@
+// This file was taken and adapted from: https://github.com/scipy/scipy/blob/main/benchmarks/asv.conf.json
+{
+    // The version of the config file format.  Do not change, unless
+    // you know what you are doing.
+    "version": 1,
+
+    // The name of the project being benchmarked
+    "project": "aeon",
+
+    // The project's homepage
+    "project_url": "https://www.aeon-toolkit.org",
+
+    // The URL of the source code repository for the project being
+    // benchmarked
+    "repo": "..",
+    "dvcs": "git",
+    "branches": ["HEAD"],
+
+    // Customizable commands for building, installing, and
+    // uninstalling the project. See asv.conf.json documentation.
+    //
+    // "install_command": ["in-dir={env_dir} python -mpip install {wheel_file}"],
+    // "uninstall_command": ["return-code=any python -mpip uninstall -y {project}"],
+    // "build_command": [
+    //     "PIP_NO_BUILD_ISOLATION=false python -m pip install . --no-deps --no-index -w {build_cache_dir} {build_dir}"
+    // ],
+
+    "build_command": [
+        "python -m build --wheel -o {build_cache_dir} {build_dir}"
+    ],
+
+    // The base URL to show a commit for the project.
+    "show_commit_url": "https://github.com/aeon-toolkit/aeon/commit",
+
+    // The Pythons you'd like to test against.  If not provided, defaults
+    // to the current version of Python used to run `asv`.
+    // "pythons": ["3.9"],
+
+    // The matrix of dependencies to test.  Each key is the name of a
+    // package (in PyPI) and the values are version numbers.  An empty
+    // list indicates to just test against the default (latest)
+    // version.
+    "matrix": {
+        "numpy": [],
+        "numba": []
+    },
+
+    // The directory (relative to the current directory) that benchmarks are
+    // stored in.  If not provided, defaults to "benchmarks"
+    "benchmark_dir": "benchmarks",
+
+    // The directory (relative to the current directory) to cache the Python
+    // environments in.  If not provided, defaults to "env"
+    "env_dir": "env",
+
+    // The tool to use to create environments.  May be "conda",
+    // "virtualenv" or other value depending on the plugins in use.
+    // If missing or the empty string, the tool will be automatically
+    // determined by looking for tools on the PATH environment
+    // variable.
+    "environment_type": "virtualenv",
+    // "environment_type": "mamba",
+    "build_cache_size": 10,
+
+    // The directory (relative to the current directory) that raw benchmark
+    // results are stored in.  If not provided, defaults to "results".
+    "results_dir": "benchmark_results",
+
+    // The directory (relative to the current directory) that the html tree
+    // should be written to.  If not provided, defaults to "html".
+    "html_dir": "html",
+
+    // The number of characters to retain in the commit hashes.
+    "hash_length": 8,
+
+    // The commits after which the regression search in `asv publish`
+    // should start looking for regressions. Dictionary whose keys are
+    // regexps matching to benchmark names, and values corresponding to
+    // the commit (exclusive) after which to start looking for
+    // regressions.  The default is to start from the first commit
+    // with results. If the commit is `null`, regression detection is
+    // skipped for the matching benchmark.
+
+    "regressions_first_commits": {
+       "io_matlab\\.StructArr\\..*": "1a002f1"
+    }
+}
diff --git a/benchmarks/benchmarks/__init__.py b/benchmarks/benchmarks/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/benchmarks/benchmarks/clustering.py b/benchmarks/benchmarks/clustering.py
new file mode 100644
index 0000000000..9032f0cc3f
--- /dev/null
+++ b/benchmarks/benchmarks/clustering.py
@@ -0,0 +1,48 @@
+from sklearn.base import BaseEstimator
+
+from .common import EstimatorBenchmark, safe_import
+
+with safe_import():
+    import aeon.clustering as aeon_clust
+
+
+class KMeansBenchmark(EstimatorBenchmark):
+    """This runs kmeans with mean averaging method."""
+
+    inits = ["random", "kmeans++"]
+
+    params = EstimatorBenchmark.params + [inits]
+    param_names = EstimatorBenchmark.param_names + ["init"]
+
+    def _build_estimator(self, k, init, distance) -> BaseEstimator:
+        return aeon_clust.TimeSeriesKMeans(
+            n_clusters=4,
+            init=init,
+            distance="euclidean",
+            averaging_method="mean",
+            n_init=1,
+            max_iter=20,
+            random_state=1,
+        )
+
+
+class KMeansBABenchmark(EstimatorBenchmark):
+    """This runs kmeans with ba averaging method."""
+
+    distances = ["dtw", "msm"]
+    params = EstimatorBenchmark.params + [
+        KMeansBenchmark.inits,
+        distances,
+    ]
+    param_names = EstimatorBenchmark.param_names + ["init", "distance"]
+
+    def _build_estimator(self, init, distance) -> BaseEstimator:
+        return aeon_clust.TimeSeriesKMeans(
+            n_clusters=4,
+            init=init,
+            distance=distance,
+            averaging_method="ba",
+            n_init=1,
+            max_iter=10,
+            random_state=1,
+        )
diff --git a/benchmarks/benchmarks/common.py b/benchmarks/benchmarks/common.py
new file mode 100644
index 0000000000..c13d693894
--- /dev/null
+++ b/benchmarks/benchmarks/common.py
@@ -0,0 +1,66 @@
+import os
+from abc import ABC, abstractmethod
+
+from sklearn.base import BaseEstimator
+
+
+class safe_import:
+
+    def __enter__(self):
+        self.error = False
+        return self
+
+    def __exit__(self, type_, value, traceback):
+        if type_ is not None:
+            self.error = True
+            suppress = not (
+                os.getenv("SCIPY_ALLOW_BENCH_IMPORT_ERRORS", "1").lower()
+                in ("0", "false")
+                or not issubclass(type_, ImportError)
+            )
+            return suppress
+
+
+class Benchmark:
+    """
+    Base class with sensible options
+    """
+
+
+with safe_import():
+    from aeon.testing.data_generation import make_example_3d_numpy
+
+
+class EstimatorBenchmark(Benchmark, ABC):
+    # Base grid (shared across all estimators)
+    shapes = [
+        (10, 1, 10),
+        (10, 1, 1000),
+        (50, 1, 100),
+        (10, 3, 10),
+        (10, 3, 1000),
+        (50, 3, 100),
+    ]
+
+    params = [shapes]
+    param_names = ["shape"]
+
+    def setup(self, shape, *est_params):
+        # Data
+        self.X_train = make_example_3d_numpy(*shape, return_y=False, random_state=1)
+        self.X_test = make_example_3d_numpy(*shape, return_y=False, random_state=2)
+
+        self.prefit_estimator = self._build_estimator(*est_params)
+        self.prefit_estimator.fit(self.X_train)
+
+    def time_fit(self, shape, *est_params):
+        est = self._build_estimator(*est_params)  # fresh each run
+        est.fit(self.X_train)
+
+    def time_predict(self, shape, *est_params):
+        self.prefit_estimator.predict(self.X_test)
+
+    @abstractmethod
+    def _build_estimator(self, *est_params) -> BaseEstimator:
+        """Return an unfitted estimator configured with the given params."""
+        ...
diff --git a/benchmarks/benchmarks/distances.py b/benchmarks/benchmarks/distances.py
new file mode 100644
index 0000000000..316126d700
--- /dev/null
+++ b/benchmarks/benchmarks/distances.py
@@ -0,0 +1,357 @@
+from abc import ABC, abstractmethod
+
+from .common import Benchmark, safe_import
+
+with safe_import():
+    import aeon.distances as aeon_dists
+    from aeon.testing.data_generation import make_example_3d_numpy
+
+
+class _DistanceBenchmark(Benchmark, ABC):
+    # params: (n_cases, n_channels, n_timepoints)
+    params = [
+        [
+            (10, 1, 10),
+            (10, 1, 1000),
+            (50, 1, 100),
+            (10, 3, 10),
+            (10, 3, 1000),
+            (50, 3, 100),
+        ],
+    ]
+    param_names = ["shape"]
+
+    def setup(self, shape):
+        self.a = make_example_3d_numpy(*shape, return_y=False, random_state=1)
+        self.b = make_example_3d_numpy(*shape, return_y=False, random_state=2)
+
+    def time_dist(self, shape):
+        self.distance_func(self.a[0], self.b[0])
+
+    def time_pairwise_dist(self, shape):
+        self.pairwise_func(self.a)
+
+    def time_one_to_multiple_dist(self, shape):
+        self.pairwise_func(self.a[0], self.b)
+
+    def time_multiple_to_multiple_dist(self, shape):
+        self.pairwise_func(self.a, self.b)
+
+    def time_muti_thread_pairwise_dist(self, shape):
+        self.pairwise_func(self.a, n_jobs=2)
+
+    def time_muti_thread_dist(self, shape):
+        self.pairwise_func(self.a, self.b, n_jobs=4)
+
+    @property
+    @abstractmethod
+    def distance_func(self):
+        """Return a callable: dist(Xi, Xj) -> float"""
+        raise NotImplementedError
+
+    @property
+    @abstractmethod
+    def pairwise_func(self):
+        """Return a callable: pairwise(X[, Y]) -> ndarray"""
+        raise NotImplementedError
+
+
+class _ElasticDistanceBenchmark(_DistanceBenchmark, ABC):
+
+    def setup(self, shape):
+        temp = make_example_3d_numpy(5, 1, 10, random_state=42, return_y=False)
+        for _ in range(3):
+            self.alignment_func(temp[0], temp[0])
+
+        super().setup(shape)
+
+    def time_alignment_path(self, shape):
+        self.alignment_func(self.a[0], self.b[0])
+
+    @property
+    @abstractmethod
+    def alignment_func(self):
+        """Return a callable: alignment_path(Xi, Xj) -> (path, cost) or similar"""
+        raise NotImplementedError
+
+
+class SquaredBenchmark(_DistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.squared_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.squared_pairwise_distance
+
+
+class ManhattanBenchmark(_DistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.manhattan_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.manhattan_pairwise_distance
+
+
+class MinkowskiBenchmark(_DistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.minkowski_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.minkowski_pairwise_distance
+
+
+class DTWBenchmark(_ElasticDistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.dtw_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.dtw_pairwise_distance
+
+    @property
+    def alignment_func(self):
+        return aeon_dists.dtw_alignment_path
+
+
+class DTWGIBenchmark(_ElasticDistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.dtw_gi_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.dtw_gi_pairwise_distance
+
+    @property
+    def alignment_func(self):
+        return aeon_dists.dtw_gi_alignment_path
+
+
+class DDTWBenchmark(_ElasticDistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.ddtw_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.ddtw_pairwise_distance
+
+    @property
+    def alignment_func(self):
+        return aeon_dists.ddtw_alignment_path
+
+
+class WDTWBenchmark(_ElasticDistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.wdtw_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.wdtw_pairwise_distance
+
+    @property
+    def alignment_func(self):
+        return aeon_dists.wdtw_alignment_path
+
+
+class WDDTWBenchmark(_ElasticDistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.wddtw_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.wddtw_pairwise_distance
+
+    @property
+    def alignment_func(self):
+        return aeon_dists.wddtw_alignment_path
+
+
+class LCSSBenchmark(_ElasticDistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.lcss_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.lcss_pairwise_distance
+
+    @property
+    def alignment_func(self):
+        return aeon_dists.lcss_alignment_path
+
+
+class ERPBenchmark(_ElasticDistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.erp_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.erp_pairwise_distance
+
+    @property
+    def alignment_func(self):
+        return aeon_dists.erp_alignment_path
+
+
+class EDRBenchmark(_ElasticDistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.edr_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.edr_pairwise_distance
+
+    @property
+    def alignment_func(self):
+        return aeon_dists.edr_alignment_path
+
+
+class TWEBenchmark(_ElasticDistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.twe_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.twe_pairwise_distance
+
+    @property
+    def alignment_func(self):
+        return aeon_dists.twe_alignment_path
+
+
+class MSMBenchmark(_ElasticDistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.msm_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.msm_pairwise_distance
+
+    @property
+    def alignment_func(self):
+        return aeon_dists.msm_alignment_path
+
+
+class ADTWBenchmark(_ElasticDistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.adtw_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.adtw_pairwise_distance
+
+    @property
+    def alignment_func(self):
+        return aeon_dists.adtw_alignment_path
+
+
+class ShapeDTWBenchmark(_ElasticDistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.shape_dtw_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.shape_dtw_pairwise_distance
+
+    @property
+    def alignment_func(self):
+        return aeon_dists.shape_dtw_alignment_path
+
+
+class SoftDTWBenchmark(_ElasticDistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.soft_dtw_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.soft_dtw_pairwise_distance
+
+    @property
+    def alignment_func(self):
+        return aeon_dists.soft_dtw_alignment_path
+
+
+class SBDBenchmark(_DistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.sbd_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.sbd_pairwise_distance
+
+
+class ShiftScaleBenchmark(_DistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.shift_scale_invariant_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.shift_scale_invariant_pairwise_distance
+
+
+class DFTSFABenchmark(_DistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.mindist_dft_sfa_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.mindist_dft_sfa_pairwise_distance
+
+
+class PAASAXBenchmark(_DistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.mindist_paa_sax_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.mindist_paa_sax_pairwise_distance
+
+
+class SAXBenchmark(_DistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.mindist_sax_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.mindist_sax_pairwise_distance
+
+
+class SFABenchmark(_DistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.mindist_sfa_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.mindist_sfa_pairwise_distance
+
+
+class MPDistBenchmark(_DistanceBenchmark):
+    @property
+    def distance_func(self):
+        return aeon_dists.mp_distance
+
+    @property
+    def pairwise_func(self):
+        return aeon_dists.mp_pairwise_distance
diff --git a/benchmarks/process_global_bencmarks.py b/benchmarks/process_global_bencmarks.py
new file mode 100644
index 0000000000..09d9fbebe4
--- /dev/null
+++ b/benchmarks/process_global_bencmarks.py
@@ -0,0 +1,49 @@
+import json
+
+import pandas as pd
+
+
+def process_global_benchmarks(f):
+    """
+    Processes the global benchmarks results into pandas DataFrame.
+
+    Parameters
+    ----------
+    f: {str, file-like}
+        Global Benchmarks output
+
+    Returns
+    -------
+    nfev, success_rate, mean_time
+        pd.DataFrame for the mean number of nfev, success_rate, mean_time
+        for each optimisation problem.
+
+    Notes
+    -----
+    Code adapted from scipy.org/scipy/benchmarks/benchmarks/process_global_benchmarks.py
+    """
+    with open(f) as fi:
+        dct = json.load(fi)
+
+    nfev = []
+    nsuccess = []
+    mean_time = []
+
+    solvers = dct[list(dct.keys())[0]].keys()
+    for _, results in dct.items():
+        _nfev = []
+        _nsuccess = []
+        _mean_time = []
+        for _, vals in results.items():
+            _nfev.append(vals["mean_nfev"])
+            _nsuccess.append(vals["nsuccess"] / vals["ntrials"] * 100)
+            _mean_time.append(vals["mean_time"])
+        nfev.append(_nfev)
+        nsuccess.append(_nsuccess)
+        mean_time.append(_mean_time)
+
+    nfev = pd.DataFrame(data=nfev, index=dct.keys(), columns=solvers)
+    nsuccess = pd.DataFrame(data=nsuccess, index=dct.keys(), columns=solvers)
+    mean_time = pd.DataFrame(data=mean_time, index=dct.keys(), columns=solvers)
+
+    return nfev, nsuccess, mean_time
diff --git a/pyproject.toml b/pyproject.toml
index b27cc4d00f..4b4cc564bf 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -97,6 +97,8 @@ dev = [
     "pytest-rerunfailures",
     "pytest-timeout",
     "pytest-xdist[psutil]",
+    "asv",
+    "snakeviz"
 ]
 binder = [
     "notebook",