diff --git a/Project.toml b/Project.toml
index f8bccd76..9f6ba6c2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ReservoirComputing"
 uuid = "7c2d2b1e-3dd4-11ea-355a-8f6a8116e294"
 authors = ["Francesco Martinuzzi"]
-version = "0.12.3"
+version = "0.12.4"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
diff --git a/README.md b/README.md
index d05a8cbc..1ee4aa68 100644
--- a/README.md
+++ b/README.md
@@ -42,6 +42,7 @@ reservoir computing models. More specifically the software offers:
     + Deep echo state networks `DeepESN`
     + Echo state networks with delayed states `DelayESN`
     + Hybrid echo state networks `HybridESN`
+    + Next generation reservoir computing `NGRC`
 - 15+ reservoir initializers and 5+ input layer initializers
 - 5+ reservoir states modification algorithms
 - Sparse matrix computation through
diff --git a/docs/pages.jl b/docs/pages.jl
index b38e1db4..b0d82b97 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -3,6 +3,7 @@ pages = [
     "Tutorials" => Any[
         "Building a model from scratch" => "tutorials/scratch.md",
         "Chaos forecasting with an ESN" => "tutorials/lorenz_basic.md",
+        "Fitting a Next Generation Reservoir Computer" => "tutorials/ngrc.md",
         #"Using Different Training Methods" => "esn_tutorials/different_training.md",
         "Deep Echo State Networks" => "tutorials/deep_esn.md",
         #"Hybrid Echo State Networks" => "tutorials/hybrid.md",
diff --git a/docs/src/api/models.md b/docs/src/api/models.md
index b75b9915..7fadf777 100644
--- a/docs/src/api/models.md
+++ b/docs/src/api/models.md
@@ -9,6 +9,16 @@
     HybridESN
 ```
 
+## Next generation reservoir computing
+
+```@docs
+    NGRC
+```
+
+```@docs
+    polynomial_monomials
+```
+
 ### Utilities
 
 ```@docs
diff --git a/docs/src/index.md b/docs/src/index.md
index dd9ae8d7..b91e3f34 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -82,6 +82,7 @@ or `dev` the package.
     + Deep echo state networks [`DeepESN`](@ref)
     + Echo state networks with delayed states [`DelayESN`](@ref)
     + Hybrid echo state networks [`HybridESN`](@ref)
+    + Next generation reservoir computing [`NGRC`](@ref)
 - 15+ reservoir initializers and 5+ input layer initializers
 - 5+ reservoir states modification algorithms
 - Sparse matrix computation through
diff --git a/docs/src/refs.bib b/docs/src/refs.bib
index 07a286dd..155f0246 100644
--- a/docs/src/refs.bib
+++ b/docs/src/refs.bib
@@ -355,3 +355,17 @@ @article{Fleddermann2025
   year = {2025},
   month = may 
 }
+
+@article{Gauthier2021,
+  title = {Next generation reservoir computing},
+  volume = {12},
+  ISSN = {2041-1723},
+  url = {http://dx.doi.org/10.1038/s41467-021-25801-2},
+  DOI = {10.1038/s41467-021-25801-2},
+  number = {1},
+  journal = {Nature Communications},
+  publisher = {Springer Science and Business Media LLC},
+  author = {Gauthier,  Daniel J. and Bollt,  Erik and Griffith,  Aaron and Barbosa,  Wendson A. S.},
+  year = {2021},
+  month = sep 
+}
diff --git a/docs/src/tutorials/ngrc.md b/docs/src/tutorials/ngrc.md
new file mode 100644
index 00000000..5e69ac2a
--- /dev/null
+++ b/docs/src/tutorials/ngrc.md
@@ -0,0 +1,171 @@
+# Next Generation Reservoir Computing
+
+This tutorial shows how to use next generation reservoir computing [NGRC](@ref)
+in ReservoirComputing.jl to model the chaotic Lorenz system.
+
+NGRC works differently compared to traditional reservoir computing. In NGRC
+the reservoir is replaced with:
+  - A delay embedding of the input
+  - A nonlinear feature map
+The model is finally trained through ridge regression, like a normal RC.
+
+
+In this tutorial we will :
+  - simulate the Lorenz system,
+  - build an NGRC model with delayed inputs and polynomial features, following the
+    [original paper](https://doi.org/10.1038/s41467-021-25801-2),
+  - train it on one-step increments,
+  - roll it out generatively and compare with the true trajectory.
+
+## 1. Setup and imports
+
+First we need to load the necessary packages. We are going to use the following:
+
+```@example ngrc
+using OrdinaryDiffEq
+using Random
+using ReservoirComputing
+using Plots
+using Statistics
+```
+
+## 2. Define Lorenz system and generate data
+
+We define the Lorenz system and integrate it to generate a long trajectory:
+
+```@example ngrc
+function lorenz!(du, u, p, t)
+    σ, ρ, β = p
+    du[1] = σ * (u[2] - u[1])
+    du[2] = u[1] * (ρ - u[3]) - u[2]
+    du[3] = u[1] * u[2] - β * u[3]
+end
+
+prob = ODEProblem(
+    lorenz!,
+    Float32[1.0, 0.0, 0.0],
+    (0.0, 200.0),
+    (10.0f0, 28.0f0, 8/3f0),
+)
+
+data = Array(solve(prob, ABM54(); dt = 0.025))  # size: (3, T)
+```
+
+We then split the time series into training and testing segments:
+
+```@example ngrc
+shift = 300
+train_len = 500
+predict_len = 900
+
+input_data = data[:, shift:(shift + train_len - 1)]
+target_data = data[:, (shift + 1):(shift + train_len)]
+test_data = data[:, (shift + train_len):(shift + train_len + predict_len - 1)]
+```
+
+## 3. Normalization
+
+It is good practice to normalize the data, especially for polynomial features:
+
+```@example ngrc
+in_mean = mean(input_data; dims = 2)
+in_std = std(input_data;  dims = 2)
+
+train_norm_x = (input_data  .- in_mean) ./ in_std
+train_norm_y = (target_data .- in_mean) ./ in_std
+test_norm_x  = (test_data   .- in_mean) ./ in_std
+
+# We train an increment (residual) model: Δy = y_{t+1} − y_t
+train_delta_y = train_norm_y .- train_norm_x
+```
+
+## 4. Build the NGRC model
+
+Now that we have the data we can start building the model.
+Following the approach of the paper we first define two feature functions:
+
+  - a constant feature
+  - a second order polynomial monomial 
+  
+```@example ngrc
+const_feature = x -> Float32[1.0]
+poly_feature  = x -> polynomial_monomials(x; degrees = 1:2)
+```
+
+Finally, we can construct the NGRC model.
+
+We set the following:
+
+```@example ngrc
+in_dims = 3
+out_dims = 3
+num_delays = 1
+```
+
+With `in_dims=3` and `num_delays=1` the delayed input length is 6.
+Adding the polinomial of degrees 1 and 2 will put give us 21 more. Finally, the constant
+term adds 1 more feature. In total we have 28 features. 
+
+We can pass the number of features to `ro_dims` to initialize the [`LinearReadout`](@ref)
+with the correct dimensions. However, unless one is planning to fry run the model without training,
+the [`train`](@ref) function will take care to adjust the dimensions.
+
+Now we build the NGRC:
+
+```@example ngrc
+rng = MersenneTwister(0)
+
+ngrc = NGRC(in_dims, out_dims; num_delays = num_delays, stride = 1, features = (const_feature, poly_feature),
+    include_input = false,  # we already encode everything in the features
+    ro_dims = 28,
+    readout_activation = identity)
+
+ps, st = setup(rng, ngrc)
+```
+
+At this point, `ngrc` is a fully specified model with:
+  - a [`DelayLayer`](@ref) that builds a 6-dimensional delayed vector from the 3D input,
+  - a [`NonlinearFeaturesLayer`](@ref) that maps that vector to 28 polynomial features,
+  - a [`LinearReadout`](@ref) (28 => 3).
+
+## 5. Training the NGRC readout
+
+We now train the linear readout using ridge regression on the increment `train_delta_y`:
+
+```@example ngrc
+ps, st = train!(ngrc, train_norm_x, train_delta_y, ps, st;
+    train_method = StandardRidge(2.5e-6))
+```
+
+where [`StandardRidge`](@ref) is the ridge regression provided natively by ReservoirComputing.jl.
+
+## 6. Generative prediction
+
+We now perform generative prediction on the increments to obtain the predicted time series:
+
+```@example ngrc
+single_step = copy(test_norm_x[:, 1]) # normalized initial condition
+traj_norm = similar(test_norm_x, 3, predict_len)
+
+for step in 1:predict_len
+    global st
+    delta_step, st = ngrc(single_step, ps, st)
+    single_step .= single_step .+ delta_step # increment update in normalized space
+    traj_norm[:, step] .= single_step
+end
+```
+
+Finally, we unscale back to the original coordinates:
+
+```@example ngrc
+traj = traj_norm .* in_std .+ in_mean # size: (3, predict_len)
+```
+
+## 7. Visualization
+
+We can now compare the predicted trajectory with the true Lorenz data on the test segment:
+
+```@example ngrc
+plot(transpose(test_data)[:, 1], transpose(test_data)[:, 2], transpose(test_data)[:, 3]; label="actual");
+plot!(transpose(traj)[:, 1], transpose(traj)[:, 2], transpose(traj)[:, 3]; label="predicted")
+```
diff --git a/src/ReservoirComputing.jl b/src/ReservoirComputing.jl
index 8c030670..ba92ea2d 100644
--- a/src/ReservoirComputing.jl
+++ b/src/ReservoirComputing.jl
@@ -40,6 +40,7 @@ include("models/esn.jl")
 include("models/esn_deep.jl")
 include("models/esn_delay.jl")
 include("models/esn_hybrid.jl")
+include("models/ngrc.jl")
 #extensions
 include("extensions/reca.jl")
 
@@ -57,9 +58,10 @@ export block_diagonal, chaotic_init, cycle_jumps, delay_line, delayline_backward
        selfloop_forwardconnection, simple_cycle, true_doublecycle
 export add_jumps!, backward_connection!, delay_line!, reverse_simple_cycle!,
        scale_radius!, self_loop!, simple_cycle!
-export train, train!, predict, resetcarry!
+export train, train!, predict, resetcarry!, polynomial_monomials
 export ESN, DeepESN, DelayESN, HybridESN
-#reca
+export NGRC
+#ext
 export RECACell, RECA
 export RandomMapping, RandomMaps
 
diff --git a/src/models/ngrc.jl b/src/models/ngrc.jl
new file mode 100644
index 00000000..0a2f8797
--- /dev/null
+++ b/src/models/ngrc.jl
@@ -0,0 +1,192 @@
+@doc raw"""
+    NGRC(in_dims, out_dims; num_delays=2, stride=1,
+         features=(), include_input=true, init_delay=zeros32,
+         readout_activation=identity, state_modifiers=(),
+         ro_dims=nothing)
+
+Next Generation Reservoir Computing (NGRC) / NVAR-style model [Gauthier2021](@cite):
+a tapped-delay embedding of the input, followed by user-defined nonlinear feature
+maps and a linear readout. This is a "reservoir-free" architecture where all dynamics
+come from explicit input delays rather than a recurrent state.
+
+`NGRC` composes:
+  1) a [`DelayLayer`](@ref) applied directly to the input, producing a vector
+     containing the current input and a fixed number of past inputs,
+  2) a [`NonlinearFeaturesLayer`](@ref) that applies user-provided functions to
+     this delayed vector and concatenates the results, and
+  3) a [`LinearReadout`](@ref) mapping the resulting feature vector to outputs.
+
+Internally, `NGRC` is represented as a [`ReservoirComputer`](@ref) with:
+  - `reservoir` = the [`DelayLayer`](@ref),
+  - `states_modifiers` = the [`NonlinearFeaturesLayer`](@ref) plus any extra
+    `state_modifiers`,
+  - `readout` = the [`LinearReadout`](@ref).
+
+## Arguments
+
+  - `in_dims`: Input dimension.
+  - `out_dims`: Output dimension.
+
+## Keyword arguments
+
+  - `num_delays`: Number of past input vectors to include. The internal
+    [`DelayLayer`](@ref) outputs a vector of length
+    `(num_delays + 1) * in_dims` (current input plus `num_delays` past inputs).
+    Default: `2`.
+  - `stride`: Delay stride in layer calls. The delay buffer is updated only when
+    the internal clock is a multiple of `stride`. Default: `1`.
+  - `init_delay`: Initializer (or tuple of initializers) for the delay history,
+    passed to [`DelayLayer`](@ref). Each initializer function is called as
+    `init(rng, in_dims, 1)` to fill one delay column. Default: `zeros32`.
+  - `features`: A function or tuple of functions `(f₁, f₂, ...)` used by
+    [`NonlinearFeaturesLayer`](@ref). Each `f` is called as `f(x)` where `x` is
+    the delayed input vector. By default it is assumed that each `f` returns a
+    vector of the same length as `x` when `ro_dims` is not provided.
+    Default: empty `()`.
+  - `include_input`: Whether to include the raw delayed input vector itself as
+    the first block of the feature vector (passed to
+    [`NonlinearFeaturesLayer`](@ref)). Default: `true`.
+  - `state_modifiers`: Extra layers applied after the `NonlinearFeaturesLayer`
+    and before the readout. Accepts a single layer, an `AbstractVector`, or a
+    `Tuple`. Default: empty `()`.
+  - `readout_activation`: Activation for the linear readout. Default: `identity`.
+  - `ro_dims`: Input dimension of the readout. If `nothing` (default), it is
+    *estimated* under the assumption that each feature function returns a
+    vector with the same length as the delayed input. In that case,
+    `ro_dims ≈ (num_delays + 1) * in_dims * n_blocks`, where `n_blocks` is the
+    number of concatenated vectors (original delayed input if
+    `include_input=true` plus one block per feature function).
+    If your feature functions change the length (e.g. constant features,
+    higher-order polynomial expansions with cross terms), you should pass
+    `ro_dims` explicitly.
+
+## Inputs
+
+  - `x :: AbstractArray (in_dims, batch)` or `(in_dims,)`
+
+## Returns
+
+  - Output `y :: (out_dims, batch)` (or `(out_dims,)` for vector input).
+  - Updated layer state (NamedTuple).
+"""
+@concrete struct NGRC <:
+                 AbstractReservoirComputer{(:reservoir, :states_modifiers, :readout)}
+    reservoir
+    states_modifiers
+    readout
+end
+
+function NGRC(in_dims::IntegerType, out_dims::IntegerType; num_delays::IntegerType = 2,
+        stride::IntegerType = 1, features = (), include_input::BoolType = True(), init_delay = zeros32,
+        readout_activation = identity, state_modifiers = (), ro_dims = nothing)
+    reservoir = DelayLayer(in_dims; num_delays = Int(num_delays), stride = Int(stride), init_delay = init_delay)
+    feats_tuple = features isa Tuple ? features : (features,)
+    nfl = NonlinearFeaturesLayer(feats_tuple...; include_input = include_input)
+    mods_tuple_raw = state_modifiers isa Tuple || state_modifiers isa AbstractVector ?
+                     (nfl, state_modifiers...) : (nfl, state_modifiers)
+    mods = _wrap_layers(mods_tuple_raw)
+    if ro_dims === nothing
+        n_taps = in_dims * (num_delays + 1)
+        inc = ReservoirComputing.known(include_input)
+        n_blocks = (inc === true ? 1 : 0) + length(feats_tuple)
+        ro_dims = n_taps * n_blocks
+        @warn """
+            NGRC: inferring readout input dimension assuming each feature f in `features`
+            returns a vector of the same length as the delayed input.
+
+            Delayed input length: n_taps = $(n_taps)
+            Blocks (input + features): $(n_blocks)
+            => ro_dims = n_taps * n_blocks = $(ro_dims)
+
+            If your feature functions change the length (e.g. constant features,
+            quadratic monomials with cross terms), please pass `ro_dims` explicitly.
+
+            Please note that, if dimensions are not correct, training will change them and
+            no error will occur.
+        """
+    end
+    readout = LinearReadout(ro_dims => out_dims, readout_activation)
+
+    return NGRC(reservoir, mods, readout)
+end
+
+function resetcarry!(
+        rng::AbstractRNG, rc::NGRC, st; init_carry = nothing)
+    carry = get(st.reservoir, :carry, nothing)
+    @warn("""
+        Next generation reservoir computing has no internal state to reset.
+        Returning untouched model states.
+        """)
+    return st
+end
+
+@doc raw"""
+    polynomial_monomials(input_vector;
+        degrees = 1:2)
+
+Generate all unordered polynomial monomials of the entries in `input_vector`
+for the given set of degrees.
+
+For each `d` in `degrees`, this function produces all degree-`d` monomials
+of the form
+
+- degree 1: `x₁, x₂, …`
+- degree 2: `x₁², x₁x₂, x₁x₃, x₂², …`
+- degree 3: `x₁³, x₁²x₂, x₁x₂x₃, x₂³, …`
+
+where combinations are taken with repetition and in non-decreasing index
+order. This means that, for example, `x₁x₂` and `x₂x₁` are represented only
+once.
+
+The returned vector is a flat list of all such products, in a deterministic
+order determined by the recursive enumeration.
+
+## Arguments
+
+- `input_vector`
+  Input vector whose entries define the variables used to build monomials.
+
+## Keyword arguments
+
+- `degrees`: An iterable of positive integers specifying which monomial degrees
+  to generate. Each degree less than `1` is skipped. Default: `1:2`.
+
+## Returns
+
+- `output_monomials` a vector of the same type as `input_vector`
+  containing all generated monomials, concatenated across the requested
+  degrees, in a deterministic order.
+"""
+function polynomial_monomials(input_vector::AbstractVector;
+        degrees = 1:2)
+    element_type = eltype(input_vector)
+    output_monomials = element_type[]
+    num_variables = length(input_vector)
+    for degree in degrees
+        degree < 1 && continue
+        index_buffer = Vector{Int}(undef, degree)
+        _polynomial_monomials_recursive!(output_monomials, input_vector,
+            index_buffer, 1, 1, num_variables
+        )
+    end
+
+    return output_monomials
+end
+
+function _polynomial_monomials_recursive!(output_monomials, input_vector,
+        index_buffer, position::Int, start_index::Int, num_variables::Int)
+    if position > length(index_buffer)
+        element_type = eltype(input_vector)
+        product_value = one(element_type)
+        @inbounds for variable_index in index_buffer
+            product_value *= input_vector[variable_index]
+        end
+        push!(output_monomials, product_value)
+    else
+        @inbounds for variable_index in start_index:num_variables
+            index_buffer[position] = variable_index
+            _polynomial_monomials_recursive!(output_monomials, input_vector,
+                index_buffer, position + 1, variable_index, num_variables)
+        end
+    end
+end
diff --git a/test/esn/test_esn.jl b/test/models/test_esn.jl
similarity index 100%
rename from test/esn/test_esn.jl
rename to test/models/test_esn.jl
diff --git a/test/esn/test_esn_deep.jl b/test/models/test_esn_deep.jl
similarity index 100%
rename from test/esn/test_esn_deep.jl
rename to test/models/test_esn_deep.jl
diff --git a/test/esn/test_esn_delay.jl b/test/models/test_esn_delay.jl
similarity index 65%
rename from test/esn/test_esn_delay.jl
rename to test/models/test_esn_delay.jl
index 9720dfbc..ffc203ee 100644
--- a/test/esn/test_esn_delay.jl
+++ b/test/models/test_esn_delay.jl
@@ -1,3 +1,9 @@
+using Test
+using Random
+using ReservoirComputing
+using Static
+using LinearAlgebra
+
 @testset "DelayESN" begin
     rng = MersenneTwister(123)
 
@@ -10,12 +16,12 @@
         desn = DelayESN(in_dims, res_dims, out_dims;
             num_delays = num_delays, stride = 1)
 
-        @test desn isa ReservoirComputer
+        @test desn isa DelayESN
 
         reservoir = desn.reservoir
         @test reservoir isa StatefulLayer
 
-        mods = desn.state_modifiers
+        mods = desn.states_modifiers
         @test mods isa Tuple
         @test !isempty(mods)
         @test first(mods) isa DelayLayer
@@ -31,28 +37,6 @@
         @test Int(ro.out_dims) == out_dims
     end
 
-    @testset "setup and forward pass shapes" begin
-        in_dims = 4
-        res_dims = 10
-        out_dims = 3
-        num_delays = 1
-
-        desn = DelayESN(in_dims, res_dims, out_dims;
-            num_delays = num_delays, stride = 1)
-
-        ps, st = setup(rng, desn)
-
-        x = rand(rng, Float32, in_dims)
-        y, st2 = desn(x, ps, st)
-        @test size(y) == (out_dims,)
-
-        X = rand(rng, Float32, in_dims, 7)
-        Y, st3 = desn(X, ps, st2)
-        @test size(Y) == (out_dims, 7)
-
-        @test propertynames(st3) == propertynames(st2)
-    end
-
     @testset "num_delays changes readout input dim" begin
         in_dims = 2
         res_dims = 6
diff --git a/test/esn/test_esn_hybrid.jl b/test/models/test_esn_hybrid.jl
similarity index 95%
rename from test/esn/test_esn_hybrid.jl
rename to test/models/test_esn_hybrid.jl
index 7f6d77c5..5aa17aa5 100644
--- a/test/esn/test_esn_hybrid.jl
+++ b/test/models/test_esn_hybrid.jl
@@ -2,6 +2,7 @@ using Test
 using Random
 using ReservoirComputing
 using Static
+using LinearAlgebra
 
 const _I32 = (m, n) -> Matrix{Float32}(I, m, n)
 const _Z32 = m -> zeros(Float32, m)
@@ -35,11 +36,11 @@ end
             state_modifiers = (),
             readout_activation = identity)
         ps, st = setup(rng, hesn)
-        @test haskey(ps, :cell) && haskey(ps, :knowledge_model) &&
+        @test haskey(ps, :reservoir) && haskey(ps, :knowledge_model) &&
               haskey(ps, :states_modifiers) && haskey(ps, :readout)
-        @test size(ps.cell.input_matrix) == (res_dims, in_dims + km_dims)
+        @test size(ps.reservoir.input_matrix) == (res_dims, in_dims + km_dims)
         @test size(ps.readout.weight) == (out_dims, res_dims + km_dims)
-        @test haskey(st, :cell) && haskey(st, :knowledge_model) &&
+        @test haskey(st, :reservoir) && haskey(st, :knowledge_model) &&
               haskey(st, :states_modifiers) && haskey(st, :readout)
     end
 
diff --git a/test/models/test_ngrc.jl b/test/models/test_ngrc.jl
new file mode 100644
index 00000000..83de3587
--- /dev/null
+++ b/test/models/test_ngrc.jl
@@ -0,0 +1,103 @@
+using Random
+using ReservoirComputing
+using LuxCore
+using Static
+
+@testset "NGRC" begin
+    rng = MersenneTwister(1234)
+
+    const_feature = x -> Float32[1.0]
+    square_feature = x -> x .^ 2
+
+    @testset "constructor & composition" begin
+        ngrc = NGRC(3, 2; num_delays = 1, stride = 2,
+            features = (const_feature, square_feature), include_input = True(),
+            init_delay = zeros32, readout_activation = tanh)
+
+        @test ngrc isa NGRC
+        @test ngrc.reservoir isa DelayLayer
+        @test ngrc.readout isa LinearReadout
+
+        dl = ngrc.reservoir
+        @test dl.in_dims == 3
+        @test dl.num_delays == 1
+        @test dl.stride == 2
+
+        @test !isempty(ngrc.states_modifiers)
+        first_mod = getfield(ngrc.states_modifiers, 1)
+        @test first_mod isa NonlinearFeaturesLayer
+    end
+
+    @testset "initialparameters & initialstates" begin
+        ngrc = NGRC(3, 2; num_delays = 1, features = (square_feature,),
+            include_input = True())
+
+        ps = initialparameters(rng, ngrc)
+        st = initialstates(rng, ngrc)
+
+        @test hasproperty(ps, :reservoir)
+        @test hasproperty(ps, :states_modifiers)
+        @test hasproperty(ps, :readout)
+
+        @test hasproperty(st, :reservoir)
+        @test hasproperty(st, :states_modifiers)
+        @test hasproperty(st, :readout)
+
+        @test ps.readout.weight isa AbstractArray
+    end
+
+    @testset "forward pass: vector input" begin
+        ngrc = NGRC(3, 2; num_delays = 1, features = (square_feature,),
+            include_input = True())
+
+        ps, st = setup(rng, ngrc)
+
+        x = rand(Float32, 3)
+        y, st2 = ngrc(x, ps, st)
+
+        @test size(y) == (2,)
+        @test propertynames(st2) == propertynames(st)
+    end
+
+    @testset "forward pass: matrix input via collectstates" begin
+        ngrc = NGRC(3, 2; num_delays = 1, features = (square_feature,),
+            include_input = True())
+
+        ps, st = setup(rng, ngrc)
+
+        X = rand(Float32, 3, 10)
+        states, st2 = collectstates(ngrc, X, ps, st)
+
+        @test size(states, 2) == size(X, 2)
+        @test size(states, 1) == ngrc.readout.in_dims
+
+        @test propertynames(st2) == propertynames(st)
+    end
+
+    @testset "simple 1D linear system learning" begin
+        # x_{t+1} = a * x_t with a = 0.8
+        a = 0.8f0
+        T = 200
+        x = zeros(Float32, T)
+        x[1] = 1.0f0
+        for t in 1:(T - 1)
+            x[t + 1] = a * x[t]
+        end
+
+        X_in = reshape(x[1:(end - 1)], 1, :)
+        Y_out = reshape(x[2:end], 1, :)
+
+        ngrc = NGRC(1, 1; num_delays = 0, stride = 1, features = (),
+            include_input = True(), ro_dims = 1)
+
+        ps, st = setup(rng, ngrc)
+
+        ps_tr, st_tr = train!(ngrc, X_in, Y_out, ps, st;
+            train_method = StandardRidge(1e-6))
+
+        @test hasproperty(ps_tr, :readout)
+        w = ps_tr.readout.weight
+        @test size(w) == (1, 1)
+        @test isapprox(w[1, 1], a; atol = 0.05)
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 5cd7aa4c..b59307c4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -13,7 +13,13 @@ end
 end
 
 @testset "Echo State Networks" begin
-    @safetestset "ESN Initializers" include("esn/test_inits.jl")
-    @safetestset "ESN model" include("esn/test_esn.jl")
-    @safetestset "DeepESN model" include("esn/test_esn_deep.jl")
+    @safetestset "ESN Initializers" include("test_inits.jl")
+    @safetestset "ESN model" include("models/test_esn.jl")
+    @safetestset "DeepESN model" include("models/test_esn_deep.jl")
+    @safetestset "DelayESN model" include("models/test_esn_delay.jl")
+    @safetestset "HybridESN model" include("models/test_esn_hybrid.jl")
+end
+
+@testset "Next Generation Reservoir Computing" begin
+    @safetestset "NGRC model" include("models/test_ngrc.jl")
 end
diff --git a/test/esn/test_inits.jl b/test/test_inits.jl
similarity index 100%
rename from test/esn/test_inits.jl
rename to test/test_inits.jl