diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 56d41cb6f82fd51f8cfd78722923f0f92a3db007..dd019b1575fd82e985cb3c876938c62b864daf8d 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -2,12 +2,7 @@
 # An example .gitlab-ci.yml file to test (and optionally report the coverage
 # results of) your [Julia][1] packages. Please refer to the [documentation][2]
 # for more information about package development in Julia.
-#
-# Here, it is assumed that your Julia package is named `MyPackage`. Change it to
-# whatever name you have given to your package.
-#
-# [1]: http://julialang.org/
-# [2]: http://julia.readthedocs.org/
+
 
 # Below is the template to run your tests in Julia
 .test_template: &test_definition
@@ -21,7 +16,7 @@
     # want coverage results.
     # - /opt/julia/bin/julia -e 'Pkg.clone(pwd()); Pkg.test("BallAndBeam",
       # coverage = true)'
-    - /opt/julia/bin/julia -e 'Pkg.checkout("BallAndBeam"); Pkg.test("BallAndBeam", coverage = true)'
+    - /opt/julia/bin/julia -e 'Pkg.test("BallAndBeam", coverage = true)'
     # Comment out below if you do not want coverage results.
     - /opt/julia/bin/julia -e 'Pkg.add("Coverage"); cd(Pkg.dir("BallAndBeam"));
       using Coverage; cl, tl = get_summary(process_folder());
diff --git a/src/BallAndBeam.jl b/src/BallAndBeam.jl
index 7c7da9c41e9fdd5e469c7287d7be21d5bce5527f..ec52790d3e654f46373403b0d3febcb717e5898c 100644
--- a/src/BallAndBeam.jl
+++ b/src/BallAndBeam.jl
@@ -18,8 +18,10 @@ export run_experiment, fra, sortfqs,
 bopl, bopl!, nypl, nypl!, plot, fbdesign, ffdesign, opendoc
 export init_sysfilter, sysfilter!, run_control_2DOF # For documentation
 
+
 using LabProcesses, Plots, Polynomials, ControlSystems, ProgressMeter
 
+include("arx.jl")
 
 
 """
diff --git a/src/arx.jl b/src/arx.jl
index 1d0569e0a5fbf32487150f2a66908625718c30c3..f670c20b1cff6404fc73c6b4410215234457986f 100644
--- a/src/arx.jl
+++ b/src/arx.jl
@@ -1,48 +1,16 @@
-abstract type Model end
-abstract type LinearModel <: Model end
-abstract type NonLinearModel <: Model end
-abstract type Network <: NonLinearModel end
-const Polynom{T<:Real} = Union{Array{T,1} , T}
-const PolynomMatrix{T} = Union{Array{Polynom{T},1},Polynom{T}, T}
+export toeplitz, getARXregressor, find_na, arx
 
 
-
-mutable struct FitStatistics
-    RMS::Float64
-    FIT::Float64
-    AIC::Float64
-    FitStatistics(y::Vector{Float64},yh::Vector{Float64},d) = new(rms(y-yh),fit(y,yh),aic(y-yh,d))
-end
-
-mutable struct FitResult
-    error::VecOrMat{Float64}
-    fit::FitStatistics
-    method::Symbol
-    λ::Float64
-    function FitResult(y,yh,d::Int,method::Symbol, λ = 0)
-        error = y-yh
-        new(error, FitStatistics(y,yh, d),method, λ)
-    end
-end
-
-mutable struct IdData
-    y::VecOrMat{Float64}
-    u::VecOrMat{Float64}
-    Ts::Float64
-    IdData(y::VecOrMat{Float64}) = new(y,[],0)
-    IdData(y::VecOrMat{Float64} ,Ts::Float64) = new(y,[],Ts)
-    IdData(y::VecOrMat{Float64}, u::VecOrMat{Float64}) = new(y,u,0)
-end
-
 ## Helper functions
 rms(x) = sqrt(mean(x.^2))
 sse(x) = sum(x.^2)
 fit(y,yh) = 100 * (1-rms(y-yh)./rms(y-mean(y)));
 aic(x,d) = log(sse(x)) + 2d/length(x)
 
-##
-Base.show(fit::FitStatistics) = println("Fit RMS:$(fit.RMS), FIT:$(fit.FIT), AIC:$(fit.AIC)")
-
+"""
+    toeplitz{T}(c::AbstractArray{T},r::AbstractArray{T})
+Returns a Toeplitz matrix where `c` is the first column and `r` is the first row.
+"""
 function toeplitz{T}(c::AbstractVector{T},r::AbstractVector{T})
     nc = length(c)
     nr = length(r)
@@ -56,19 +24,25 @@ function toeplitz{T}(c::AbstractVector{T},r::AbstractVector{T})
 end
 
 """
-`a::Vector{Float64}`: The polynomial coeffs A(z) (not including the first 1)\n
-`b::VecOrMat{Float64}`: The polynomial coeffs B(z)\n
-`na::Int`\n
-`nb::Vector{Int}`\n
+    getARXregressor(y::AbstractVector,u::AbstractVecOrMat, na, nb)
+Returns a shortened output signal `y` and a regressor matrix `A` such that the least-squares ARX model estimate of order `na,nb` is `y\\A`
+Return a regressor matrix used to fit an ARX model on, e.g., the form
+`A(z)y = B(z)f(u)`
+with output `y` and input `u` where the order of autoregression is `na` and
+the order of input moving average is `nb`
+# Example
+Here we test the model with the Function `f(u) = √(|u|)`
+```julia
+A     = [1,2*0.7*1,1] # A(z) coeffs
+B     = [10,5] # B(z) coeffs
+u     = randn(100) # Simulate 100 time steps with Gaussian input
+y     = filt(B,A,u)
+yr,A  = getARXregressor(y,u,3,2) # We assume that we know the system order 3,2
+x     = A\yr # Estimate model polynomials
+plot([yr A*x], lab=["Signal" "Prediction"])
+```
+For nonlinear ARX-models, see [BasisFunctionExpansions.jl](https://github.com/baggepinnen/BasisFunctionExpansions.jl/)
 """
-type ARX <: LinearModel
-    a::Polynom
-    b::PolynomMatrix
-    na::Int
-    nb::Polynom{Int}
-end
-
-
 function getARXregressor(y::AbstractVector{Float64},u::AbstractVecOrMat{Float64}, na, nb)
     assert(length(nb) == size(u,2))
     m    = max(na+1,maximum(nb))
@@ -84,13 +58,11 @@ function getARXregressor(y::AbstractVector{Float64},u::AbstractVecOrMat{Float64}
     return y,A
 end
 
-getARXregressor(iddata::IdData, na, nb) = getARXregressor(iddata.y,iddata.u, na, nb)
-
 """
-    find_na(y::Vector,n::Int)
+    find_na(y::AbstractVector,n::Int)
 Plots the RMSE and AIC for model orders up to `n`. Useful for model selection
 """
-function find_na(y::Vector,n::Int)
+function find_na(y::AbstractVector,n::Int)
     error = zeros(n,2)
     for i = 1:n
         w,e = ar(y,i)