diff --git a/src/SystemIdentification.jl b/src/SystemIdentification.jl
index 6282c6e0afd2a88e6b8c05bd7b8c92c5483d818a..0264d4f597e0854938e11e08858704ff4a323672 100644
--- a/src/SystemIdentification.jl
+++ b/src/SystemIdentification.jl
@@ -11,10 +11,12 @@ FitStatistics,
 FitResult,
 IdData,
 # Functions
-ar,arx,getARregressor,getARXregressor,find_na
+ar,arx,getARregressor,getARXregressor,find_na,
+toeplitz, kalman
 
 ## Fit Methods =================
 :LS
+:LS_reg
 :LM
 
 ## Types =======================
@@ -22,8 +24,8 @@ abstract Model
 abstract LinearModel <: Model
 abstract NonLinearModel <: Model
 abstract Network <: NonLinearModel
-typealias Polynom Union(Real, Array{Real,1})
-typealias PolynomMatrix Union(Array{Polynom,1},Polynom)
+typealias Polynom{T<:Real} Union(T, Array{T,1})
+typealias PolynomMatrix{T<:Real} Union(Array{Polynom,1},Polynom)
 
 
 """
@@ -32,7 +34,7 @@ typealias PolynomMatrix Union(Array{Polynom,1},Polynom)
 `bias::Bool`\n
 `λ::Float64`\n
 """
-type AR <: Model
+type AR <: LinearModel
     a::Polynom
     na::Int
     bias::Bool
@@ -47,7 +49,7 @@ end
 `bias::Bool`\n
 `λ::Float64`\n
 """
-type ARX <: Model
+type ARX <: LinearModel
     a::Polynom
     b::PolynomMatrix
     na::Int
@@ -68,14 +70,17 @@ type FitStatistics
     RMS::Float64
     FIT::Float64
     AIC::Float64
-    FitStatistics(e::Vector{Float64}) = new(rms(e),fit(e),aic(e))
+    FitStatistics(y::Vector{Float64},yh::Vector{Float64},d) = new(rms(y-yh),fit(y,yh),aic(y-yh,d))
 end
 
 type FitResult
     error::VecOrMat{Float64}
     fit::FitStatistics
     method::Symbol
-    FitResult(error::VecOrMat{Float64},method::Symbol) = new(error, FitStatistics(error),method)
+    function FitResult(y,yh,d::Int,method::Symbol)
+        error = y-yh
+        new(error, FitStatistics(y,yh, d),method)
+    end
 end
 
 type IdData
@@ -87,7 +92,11 @@ type IdData
     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)
 
 include("armax.jl")
 include("kalman.jl")
diff --git a/src/armax.jl b/src/armax.jl
index b3aafce6ab015298fb70bebe8448742c53beb5d4..7f1e97b519732bfda440dfe0ab0975a11f3f91dc 100644
--- a/src/armax.jl
+++ b/src/armax.jl
@@ -9,7 +9,10 @@ function ar(y::Vector{Float64}, na; λ = 0, bias=false)
     end
     prediction = A*w
     error = y_train - prediction
-    return AR(w,na,bias,λ), FitResult(error,:LS)
+    model = AR(w,na,bias,λ)
+    na = bias ? na+1:na
+    result = FitResult(y_train,prediction,na, λ>0?(:LS_reg) :(:LS))
+    return model, result
 end
 ar(iddata::IdData, na; λ = 0, doplot=false, bias=false) = ar(iddata.y, na; λ = 0, bias=bias)
 
@@ -24,8 +27,10 @@ function arx(y::Vector{Float64}, u::VecOrMat{Float64}, na, nb; λ = 0, doplot=fa
     end
     prediction = A*w
     error = y_train - prediction
-
-    return AR(w,na,bias,λ), FitResult(error,:LS)
+    model = AR(w,na,bias,λ)
+    na = bias ? na+1:na
+    result = FitResult(y_train,prediction,na, λ>0?(:LS_reg) :(:LS))
+    return model, result
 end
 arx(iddata::IdData, na, nb; λ = 0,  bias=false) = arx(iddata.y, iddata.u, na, nb; λ = 0, bias=bias)
 
diff --git a/test/tests.jl b/test/tests.jl
index d8dd67691b864522056f065967921bbed068ffe4..4ded4dbf49dfccd6e7746926f24d7e5064124731 100644
--- a/test/tests.jl
+++ b/test/tests.jl
@@ -1,5 +1,8 @@
+using SystemIdentification
 @assert isa(1,Polynom)
 @assert isa(1.0,Polynom)
 @assert isa([1.0; 1.0],Polynom)
 @assert isa(1.0,PolynomMatrix)
 @assert isa(1,PolynomMatrix)
+
+ar(collect(1:5.0),1,bias=true)