diff --git a/Project.toml b/Project.toml
index 77c55a4..48269b4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "SolveDSGE"
 uuid = "00bf1f32-23ad-54cc-bf6e-3216db8a43a2"
 authors = ["Richard Dennis <richard.dennis@glasgow.ac.uk>"]
-version = "0.4.21"
+version = "0.4.22"
 
 [deps]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
diff --git a/src/analysis_functions.jl b/src/analysis_functions.jl
index 8a65bce..003f0b4 100644
--- a/src/analysis_functions.jl
+++ b/src/analysis_functions.jl
@@ -2943,6 +2943,10 @@ end
 
 function euler_errors(model::REModel,soln::R,domain::Union{Array{T,2},Array{T,1}},npoints::S,seed::S = 123456) where {S<:Integer,T<:AbstractFloat,R<:PerturbationSolutionDet}
 
+    if model.solvers == "Perturbation"
+        error("The model hase been restricted to perturbation solvers only")
+    end
+
     nx = length(soln.hbar)
     ny = length(soln.gbar)
     f = zeros(nx+ny)
@@ -2966,6 +2970,10 @@ end
 
 function euler_errors(model::REModel,soln::R,domain::Union{Array{T,2},Array{T,1}},npoints::S,seed::S = 123456) where {S<:Integer,T<:AbstractFloat,R<:PerturbationSolutionStoch}
 
+    if model.solvers == "Perturbation"
+        error("The model hase been restricted to perturbation solvers only")
+    end
+
     nx = length(soln.hbar)
     ny = length(soln.gbar)
     ns = size(soln.k,2)
@@ -2992,6 +3000,10 @@ end
 
 function euler_errors(model::REModel,soln::R,npoints::S,seed::S = 123456) where {S<:Integer,R<:ProjectionSolutionDet}
 
+    if model.solvers == "Perturbation"
+        error("The model hase been restricted to perturbation solvers only")
+    end
+
     nx = size(soln.domain,2)
     nv = length(soln.variables)
     ny = nv - nx
@@ -3017,6 +3029,10 @@ end
 
 function euler_errors(model::REModel,soln::R,npoints::S,seed::S = 123456) where {S<:Integer,R<:ProjectionSolutionStoch}
 
+    if model.solvers == "Perturbation"
+        error("The model hase been restricted to perturbation solvers only")
+    end
+
     nx = size(soln.domain,2)
     nv = length(soln.variables)
     ny = nv - nx
@@ -3122,21 +3138,21 @@ function check_taylor_convergence(model::REModel,ss::Array{T,1},z::Array{T,1},de
             ddddddd(x) = ForwardDiff.jacobian(dddddd,x,ForwardDiff.JacobianConfig(dddddd,x,ForwardDiff.Chunk{2}()))[:,1:2*nv]
 
             for j = 1:degree
-                kprod = kron(kprod,(ss - z))
+                kprod = kron(kprod,(z - ss))
                 if j == 1
-                    taylor_terms[i,j+1] = d(point)'kprod
+                    taylor_terms[i,j+1] = abs.(d(point))'kprod
                 elseif j == 2
-                    taylor_terms[i,j+1] = vec(dd(point))'kprod/factorial(j)
+                    taylor_terms[i,j+1] = abs.(vec(dd(point)))'kprod/factorial(j)
                 elseif j == 3
-                    taylor_terms[i,j+1] = vec(ddd(point))'kprod/factorial(j)
+                    taylor_terms[i,j+1] = abs.(vec(ddd(point)))'kprod/factorial(j)
                 elseif j == 4
-                    taylor_terms[i,j+1] = vec(dddd(point))'kprod/factorial(j)
+                    taylor_terms[i,j+1] = abs.(vec(dddd(point)))'kprod/factorial(j)
                 elseif j == 5
-                    taylor_terms[i,j+1] = vec(ddddd(point))'kprod/factorial(j)
+                    taylor_terms[i,j+1] = abs.(vec(ddddd(point)))'kprod/factorial(j)
                 elseif j ==6
-                    taylor_terms[i,j+1] = vec(dddddd(point))'kprod/factorial(j)
+                    taylor_terms[i,j+1] = abs.(vec(dddddd(point)))'kprod/factorial(j)
                 else
-                    taylor_terms[i,j+1] = vec(ddddddd(point))'kprod/factorial(j)
+                    taylor_terms[i,j+1] = abs.(vec(ddddddd(point)))'kprod/factorial(j)
                 end
             end
         end
@@ -3183,24 +3199,24 @@ function check_taylor_convergence(model::REModel,ss::Array{T,1},z::Array{T,1},de
             ddddddd(x) = ForwardDiff.jacobian(dddddd,x,ForwardDiff.JacobianConfig(dddddd,x,ForwardDiff.Chunk{2}()))[:,1:2*nv]
 
             for j = 1:degree
-                kprod = kron(kprod,(ss - z))
+                kprod = kron(kprod,(z - ss))
                 if j == 1
-                    taylor_terms[i,j+1] = d(point)'kprod
+                    taylor_terms[i,j+1] = abs.(d(point))'kprod
                 elseif j == 2
-                    taylor_terms[i,j+1] = vec(dd(point))'kprod/factorial(j)
+                    taylor_terms[i,j+1] = abs.(vec(dd(point)))'kprod/factorial(j)
                 elseif j == 3
-                    taylor_terms[i,j+1] = vec(ddd(point))'kprod/factorial(j)
+                    taylor_terms[i,j+1] = abs.(vec(ddd(point)))'kprod/factorial(j)
                 elseif j == 4
-                    taylor_terms[i,j+1] = vec(dddd(point))'kprod/factorial(j)
+                    taylor_terms[i,j+1] = abs.(vec(dddd(point)))'kprod/factorial(j)
                 elseif j == 5
-                    taylor_terms[i,j+1] = vec(ddddd(point))'kprod/factorial(j)
+                    taylor_terms[i,j+1] = abs.(vec(ddddd(point)))'kprod/factorial(j)
                 elseif j ==6
-                    taylor_terms[i,j+1] = vec(dddddd(point))'kprod/factorial(j)
+                    taylor_terms[i,j+1] = abs.(vec(dddddd(point)))'kprod/factorial(j)
                 else
-                    taylor_terms[i,j+1] = vec(ddddddd(point))'kprod/factorial(j)
+                    taylor_terms[i,j+1] = abs.(vec(ddddddd(point)))'kprod/factorial(j)
                 end
             end
         end
         return taylor_terms
     end
-end
\ No newline at end of file
+end
diff --git a/src/parser_functions.jl b/src/parser_functions.jl
index 22e19d0..67f6176 100644
--- a/src/parser_functions.jl
+++ b/src/parser_functions.jl
@@ -641,6 +641,17 @@ function create_steady_state_equations(model::ModelPrimatives)
 
     sorted_combined_names = combined_names[sortperm(length.(combined_names),rev = true)]
 
+        #= First we go through every equation and replace exp with : and log with ;.  
+       This is to guard them during variables and parameter substitution. =#
+
+       for i = 1:length(steady_state_equations)
+        if occursin("exp",equations[i]) == true
+            steady_state_equations[i] = replace(steady_state_equations[i],"exp" => ":")
+        elseif occursin("log",steady_state_equations[i]) == true
+            steady_state_equations[i] = replace(steady_state_equations[i],"log" => ";")
+        end
+    end
+
     # Now we go through every equation and replace future variables, variables, and
     # shocks with a numbered element of a vector, "x".  We also replace parameter
     # names with parameter values
@@ -690,6 +701,16 @@ function create_steady_state_equations(model::ModelPrimatives)
         end
     end
 
+    #= Finally, go back through every equation and restore exp and log where necessary =#
+
+    for i = 1:length(steady_state_equations)
+        if occursin(":",steady_state_equations[i]) == true
+            steady_state_equations[i] = replace(steady_state_equations[i],":" => "exp")
+        elseif occursin(";",steady_state_equations[i]) == true
+            steady_state_equations[i] = replace(steady_state_equations[i],";" => "log")
+        end
+    end
+
     return steady_state_equations
 
 end
diff --git a/src/structures.jl b/src/structures.jl
index df65eb8..3ae1740 100644
--- a/src/structures.jl
+++ b/src/structures.jl
@@ -5,7 +5,7 @@ abstract type ModelPrimatives end
 
 struct REModelPrimatives{Q<:AbstractString} <: ModelPrimatives
 
-    # This structure contains the critial model information extracted from a
+    # This structure contains the crucial model information extracted from a
     # model file.
 
     states::Array{Q,1}
@@ -497,4 +497,4 @@ struct DenHaanMarcetStatistic{T<:Real,S<:Integer}
     ten_percent::T
     degrees_of_freedom::S
 
-end
\ No newline at end of file
+end