diff --git a/deps/build.jl b/deps/build.jl
index 481c66d1bd9fe9f515ca8bccf89bebd0e42c15a7..8838def3385d2997bc04d71a8753ec98eeb27c12 100644
--- a/deps/build.jl
+++ b/deps/build.jl
@@ -3,7 +3,7 @@ using BinaryProvider
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
-PREPROCESSOR_VERSION = "46a2272ef929a7acbf276fea2c43089bfc7d0e34"
+PREPROCESSOR_VERSION = "70ec52a6cce2ad116b0de68eec505a4092e50a82"
REMOTE_PATH = "https://dynare.adjemian.eu/preprocessor/$PREPROCESSOR_VERSION"
products = Product[
@@ -11,11 +11,11 @@ products = Product[
]
download_info = Dict(
- Linux(:i686, :glibc) => ("$REMOTE_PATH/linux/32/preprocessor.tar.gz", "3e847abfd259e3d21a28a315da322f3b326910e8f1d198efdc0bd200eaea00ec"),
- Linux(:x86_64, :glibc) => ("$REMOTE_PATH/linux/64/preprocessor.tar.gz", "bd31db8a0b29f43ba63d8df7a96cde5349851bd72c562200dbfa9d37f5606219"),
- MacOS() => ("$REMOTE_PATH/osx/64/preprocessor.tar.gz", "001fa0ab957623b1e5b01b8b7563b92f960adda79027f054d1f72df8ad60c30b"),
- Windows(:i686) => ("$REMOTE_PATH/windows/32/preprocessor.tar.gz", "1be8d8fc8cc798dec811f36525ff5b5fd6bf78e84be659068df30d4cf737eaac"),
- Windows(:x86_64) => ("$REMOTE_PATH/windows/64/preprocessor.tar.gz", "32485a1dde052e97b7844d731af1f1c0c47c5154f22e29b71789ffe9a1f25e8c"),
+ Linux(:i686, :glibc) => ("$REMOTE_PATH/linux/32/preprocessor.tar.gz", "dcfc2a15c35a6ee916cd0b20cec89a454082ad3aa512ea2dfd62ea2a22e58561"),
+ Linux(:x86_64, :glibc) => ("$REMOTE_PATH/linux/64/preprocessor.tar.gz", "1f9c619c2116542c0810ae2ba8d12525993d03a46d53322e49b3fad35f6c51e7"),
+ MacOS() => ("$REMOTE_PATH/osx/64/preprocessor.tar.gz", "4d082dbbee46810b7350d9594f1c1d848e37c871f2fa7456be60f2847682d3bf"),
+ Windows(:i686) => ("$REMOTE_PATH/windows/32/preprocessor.tar.gz", "e277c9dc9248ff57476612e42d7fd691ffd72cee1de4aea4926544906cc772e3"),
+ Windows(:x86_64) => ("$REMOTE_PATH/windows/64/preprocessor.tar.gz", "5a8cbe1f58664446c28e074b2e5b79aa23c64ccd84c7b0425d2003a29292af64"),
)
for p in products
diff --git a/src/DynareSolvers.jl b/src/DynareSolvers.jl
new file mode 100644
index 0000000000000000000000000000000000000000..031e5390873b9cee538ee4e3954ff75c8700283f
--- /dev/null
+++ b/src/DynareSolvers.jl
@@ -0,0 +1,248 @@
+module DynareSolvers
+
+##
+ # Copyright (C) 2018 Dynare Team
+ #
+ # This file is part of Dynare.
+ #
+ # Dynare is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU General Public License as published by
+ # the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # Dynare is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU General Public License
+ # along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+##
+
+export trustregion
+
+"""
+ dogleg!(x::Vector{Float64}, r::Matrix{Float64}, b::Vector{Float64}, d::Vector{Float64}, δ::Float64)
+
+Given an `n` by `n` matrix `r`, an `n` by 1 vector `d` with non zero entries, an `n` by `1`
+vector `b`, and a positive number δ, the problem is to determine the convex combination `x`
+of the gauss-newton and scaled gradient directions that minimizes (r*x - b) in the least
+squares sense, subject to the restriction that the euclidean norm of d*x be at most delta.
+"""
+function dogleg!(x::Vector{Float64}, r::Matrix{Float64}, b::Vector{Float64}, d::Vector{Float64}, δ::Float64)
+ n = length(x)
+ @assert length(d)==n
+ @assert length(b)==n
+ @assert size(r,1)==n
+ @assert size(r,2)==n
+ # Compute the Gauss-Newton direction.
+ x .= r\b
+ # Compute norm of scaled x.
+ qnorm = zero(Float64)
+ @inbounds for i = 1:n
+ qnorm += (d[i]*x[i])^2
+ end
+ qnorm = sqrt(qnorm)
+ if qnorm<=δ
+ # Gauss-Newton direction is acceptable. There is nothing to do here.
+ else
+ # Gauss-Newton direction is not acceptable…
+ # Compute the scale gradient direction.
+ s = (r'*b)./d
+ # Compute the norm of the scaled gradient direction.
+ gnorm = norm(s)
+ if gnorm>0
+ # Normalize and rescale → gradient direction.
+ @inbounds for i = 1:n
+ s[i] = s[i]/(gnorm*d[i])
+ end
+ temp = norm(r*s)
+ sgnorm = gnorm/(temp*temp)
+ if sgnorm<=δ
+ # The scaled gradient direction is not acceptable…
+ # Compute the point along the dogleg at which the
+ # quadratic is minimized.
+ bnorm = norm(b)
+ temp1 = δ/qnorm
+ temp2 = sgnorm/δ
+ temp0 = bnorm*bnorm*temp2/(gnorm*qnorm)
+ temp0 -= temp1*temp2*temp2-sqrt((temp0-temp1)^2+(one(Float64)-temp1*temp1)*(one(Float64)-temp2*temp2))
+ α = temp1*(one(Float64)-temp2*temp2)/temp0
+ else
+ # The scaled gradient direction is acceptable.
+ α = zero(Float64)
+ end
+ else
+ # If the norm of the scaled gradient direction is zero.
+ α = δ/qnorm
+ sgnorm = .0
+ end
+ # Form the appropriate convex combination of the Gauss-Newton direction and the
+ # scaled gradient direction.
+ x .= α*x + (one(Float64)-α)*min(sgnorm, δ)*s
+ end
+end
+
+"""
+ trustregion(fj::Function, y0::Vector{Float64}, gstep::Float64, tolf::Float64, tolx::Float64, maxiter::Int)
+
+Solves a system of nonlinear equations of several variables using a trust region method.
+"""
+function trustregion(f!::Function, j!::Function, x0::Vector{Float64}, tolx::Float64, tolf::Float64, maxiter::Int)
+ macheps = eps(Float64)
+ n, iter, info = length(x0), 1, 0
+ p1, p5, p001, p0001 = .1, .5, .001, .0001
+ t1, t2, t3, t4 = .1, .5, .001, 1.0e-4
+ fnorm, fnorm1 = one(Float64), one(Float64)
+ x = copy(x0)
+ xx = Vector{Float64}(n)
+ fval0 = Vector{Float64}(n) # residuals of the non linear equations
+ fval1 = Vector{Float64}(n) # residuals of the non linear equations (next)
+ fjac = Matrix{Float64}(n,n) # jacobian matrix of the non linear equations
+ fjaccnorm = Vector{Float64}(n) # norms of the columns of the Jacobian matrix
+ fjaccnorm__ = Vector{Float64}(n) # copy of fjaccnorm
+ w = Vector{Float64}(n)
+ p = Vector{Float64}(n)
+ # Initial evaluation of the residuals (and compute the norm of the residuals)
+ try
+ f!(fval0, x)
+ fnorm = norm(fval0)
+ catch
+ error("The system of nonlinear equations cannot be evaluated on the initial guess!")
+ end
+ # Initial evaluation of the Jacobian
+ try
+ j!(fjac, x)
+ catch
+ error("The Jacobian of the system of nonlinear equations cannot be evaluated on the initial guess!")
+ end
+ # Initialize counters
+ ncsucc = zero(Int)
+ nslow1 = zero(Int)
+ nslow2 = zero(Int)
+ # Newton iterations
+ while iter<=maxiter && info==0
+ # Compute columns norm for the Jacobian matrix.
+ @inbounds for i=1:n
+ fjaccnorm[i] = zero(Float64)
+ @inbounds for j = 1:n
+ fjaccnorm[i] += fjac[j,i]*fjac[j,i]
+ end
+ fjaccnorm[i] = sqrt(fjaccnorm[i])
+ end
+ if iter==1
+ # On the first iteration, calculate the norm of the scaled vector of unknwonws x
+ # and initialize the step bound δ. Scaling is done according to the norms of the
+ # columns of the initial jacobian.
+ @inbounds for i = 1:n
+ fjaccnorm__[i] = fjaccnorm[i]
+ end
+ fjaccnorm__[find(abs.(fjaccnorm__).<macheps)] = one(Float64)
+ xnorm = zero(Float64)
+ @inbounds for i=1:n
+ xnorm += (fjaccnorm__[i]*x[i])^2
+ end
+ xnorm = sqrt(xnorm)
+ δ = max(xnorm, one(Float64))
+ else
+ fjaccnorm__ = max.(.1*fjaccnorm__, fjaccnorm)
+ end
+ # Determine the direction p (with trust region model defined in dogleg routine).
+ dogleg!(p, fjac, fval0, fjaccnorm__, δ)
+ @inbounds for i=1:n
+ w[i] = fval0[i]
+ @inbounds for j = 1:n
+ w[i] -= fjac[i,j]*p[j]
+ end
+ end
+ # Compute the norm of p.
+ pnorm = zero(Float64)
+ @inbounds for i=1:n
+ pnorm += (fjaccnorm__[i]*p[i])^2
+ end
+ pnorm = sqrt(pnorm)
+ # On first iteration adjust the initial step bound.
+ if iter==1
+ δ = min(δ, pnorm)
+ end
+ # Evaluate the function at xx = x+p and calculate its norm.
+ @inbounds for i = 1:n
+ xx[i] = x[i]-p[i]
+ end
+ f!(fval1, xx)
+ fnorm1 = norm(fval1)
+ # Compute the scaled actual reduction.
+ if fnorm1<fnorm
+ actualreduction = one(Float64)-(fnorm1/fnorm)^2
+ else
+ actualreduction = -one(Float64)
+ end
+ # Compute the scaled predicted reduction and the ratio of the actual to the
+ # predicted reduction.
+ t = norm(w)
+ ratio = zero(Float64)
+ if t<fnorm
+ predictedreduction = one(Float64) - (t/fnorm)^2
+ ratio = actualreduction/predictedreduction
+ end
+ # Update the radius of the trust region if need be.
+ if ratio<p1
+ # Reduction is much smaller than predicted… Reduce the radius of the trust region.
+ ncsucc = 0
+ δ = p5*δ
+ else
+ ncsucc += 1
+ if ratio>=p5 || ncsucc>1
+ δ = max(δ,pnorm/p5)
+ end
+ if abs(ratio-one(Float64))<p1
+ δ = pnorm/p5
+ end
+ end
+ if ratio>=1.0e-4
+ # Succesfull iteration. Update x, xnorm, fval0, fnorm and fjac.
+ xnorm = zero(Float64)
+ @inbounds for i=1:n
+ x[i] = xx[i]
+ xnorm += (fjaccnorm__[i]*x[i])^2
+ fval0[i] = fval1[i]
+ end
+ xnorm = sqrt(xnorm)
+ fnorm = fnorm1
+ end
+ # Determine the progress of the iteration.
+ nslow1 += 1
+ if actualreduction>=p001
+ nslow1 = 0
+ end
+ # Test for convergence.
+ if δ<tolx*xnorm || fnorm<tolf
+ info = 1
+ continue
+ end
+ # Tests for termination and stringent tolerances.
+ if p1*max(p1*δ, pnorm)<=macheps*xnorm
+ # xtol is too small. no further improvement in
+ # the approximate solution x is possible.
+ info = 3
+ continue
+ end
+ if nslow1==15
+ # iteration is not making good progress, as
+ # measured by the improvement from the last
+ # fifteen iterations.
+ info = 4
+ continue
+ end
+ # Compute the jacobian for next the iteration.
+ j!(fjac, x)
+ iter += 1
+ end
+ if info==0 && iter>maxiter
+ info = 2
+ fill!(x, Inf)
+ end
+ return x, info
+end
+
+end
diff --git a/src/SteadyState.jl b/src/SteadyState.jl
index ab68c52955b2e7ed4dfbc609eac6b57fb749b346..2c4fbc7a46df83956b58c32ac720b8cdb475af66 100644
--- a/src/SteadyState.jl
+++ b/src/SteadyState.jl
@@ -1,7 +1,7 @@
module SteadyState
##
- # Copyright (C) 2016 Dynare Team
+ # Copyright (C) 2016-2018 Dynare Team
#
# This file is part of Dynare.
#
@@ -19,17 +19,19 @@ module SteadyState
# along with Dynare. If not, see <http://www.gnu.org/licenses/>.
##
-using NLsolve
+using DynareSolvers
+using AutoAligns
import DynareModel.Model
import DynareOutput.Output
+export staticmodule
export steady, steady!
export steady_state, steady_state!
function steady(model::Model, oo::Output)
if model.analytical_steady_state || model.user_written_analytical_steady_state
- steadystate = zeros(length(model.endo))
+ steadystate = zeros(model.endo_nbr)
model.steady_state(steadystate, oo.exo_steady_state, model.params)
return steadystate
else
@@ -40,61 +42,109 @@ end
function steady!(model::Model, oo::Output)
if model.analytical_steady_state || model.user_written_analytical_steady_state
model.steady_state(oo.steady_state, oo.exo_steady_state, model.params)
- return
else
error("You have to provide a closed form solution for the steady state, or declare a guess\nfor the steady state as a third input argument.")
end
end
function steady(model::Model, oo::Output, yinit::Vector{Float64})
- f!(fval::Vector{Float64}, y::Vector{Float64}) = model.static(y, oo.exo_steady_state, model.params, fval)
- j!(fjac::Matrix{Float64}, y::Vector{Float64}) = model.static(y, oo.exo_steady_state, model.params, fjac)
- fj!(fval::Vector{Float64}, fjac::Matrix{Float64}, y::Vector{Float64}) = model.static(y, oo.exo_steady_state, model.params, fval, fjac)
- r = nlsolve(f!, j!, fj!, yinit, show_trace=false)
- if converged(r)
- return r.zero
+ T = zeros(Float64, sum(model.temporaries.static[1:2]))
+ f!(fval::Vector{Float64}, y::Vector{Float64}) = model.static(T, fval, y , oo.exo_steady_state, model.params)
+ j!(fjac::Matrix{Float64}, y::Vector{Float64}) = model.static(T, fjac, y, oo.exo_steady_state, model.params)
+ ys, info = trustregion(f!, j!, yinit, 1.0e-6, 1.0e-6, 100)
+ return ys, info
+end
+
+function steady!(model::Model, oo::Output, yinit::Vector{Float64})
+ T = zeros(Float64, sum(model.temporaries.static[1:2]))
+ f!(fval::Vector{Float64}, y::Vector{Float64}) = model.static(T, fval, y , oo.exo_steady_state, model.params)
+ j!(fjac::Matrix{Float64}, y::Vector{Float64}) = model.static(T, fjac, y, oo.exo_steady_state, model.params)
+ oo.steady_state, info = trustregion(f!, j!, yinit, 1.0e-6, 1.0e-6, 100)
+end
+
+function steady_state(model::Model, oo::Output, yinit::Vector{Float64}, display::Bool=true)
+ ys, info = steady(model, oo, yinit)
+ if info==1
+ if display
+ display_steady_state(model, oo, ys)
+ end
else
- return fill(NaN, length(yinit))
+ error("Steady state not found!")
end
+ return ys, info
end
-function steady!(model::Model, oo::Output, yinit::Vector{Float64})
- f!(fval::Vector{Float64}, y::Vector{Float64}) = model.static(y, oo.exo_steady_state, model.params, fval)
- j!(fjac::Matrix{Float64}, y::Vector{Float64}) = model.static(y, oo.exo_steady_state, model.params, fjac)
- fj!(fval::Vector{Float64}, fjac::Matrix{Float64}, y::Vector{Float64}) = model.static(y, oo.exo_steady_state, model.params, fval, fjac)
- r = nlsolve(f!, j!, fj!, yinit, show_trace=false)
- if converged(r)
- oo.steady_state = r.zero
+function steady_state!(model::Model, oo::Output, yinit::Vector{Float64}, display::Bool=true)
+ steady!(model, oo, yinit)
+ if issteadystate(model, oo, oo.steady_state)
+ if display
+ display_steady_state(model, oo)
+ end
else
- oo.steady_state = fill(NaN, length(yinit))
+ error("Steady state not found!")
end
end
-function steady_state(model::Model, oo::Output)
+function steady_state(model::Model, oo::Output, display::Bool=true)
ys = steady(model, oo)
- display_steady_state(model, oo, ys)
+ if !issteadystate(model, oo, ys)
+ error("Steady state provided in steady state block/file is not correct!")
+ end
+ if display
+ display_steady_state(model, oo, ys)
+ end
+ return ys, info
end
-function steady_state!(model::Model, oo::Output)
+function steady_state!(model::Model, oo::Output, display::Bool=true)
steady!(model, oo)
- display_steady_state(model, oo, oo.steady_state)
+ if !issteadystate(model, oo, oo.steady_state)
+ error("Steady state provided in steady state block/file is not correct!")
+ end
+ if issteadystate(oo.steady_state)
+ if display
+ display_steady_state(model, oo)
+ end
+ else
+ error("Steady state not found!")
+ end
end
-function display_steady_state(model::Model, oo::Output, ys::Vector{Float64})
- println("\n\nSTEADY STATE:\n")
- for i = 1:length(model.endo)
- println(string(model.endo[i].name, " = ", ys[i]))
+function display_steady_state(model::Model, oo::Output)
+ aa = AutoAlign(align = Dict(1 => left, 2=> left, :default => right))
+ println("\n Deterministic Steady State:\n")
+ for i = 1:model.endo_nbr
+ if abs(oo.steady_state[i])<1.0e-6
+ tmp = zero(Float64)
+ else
+ tmp = oo.steady_state[i]
+ end
+ print(aa, " ", model.endo[i].name, "\t=\t", @sprintf("%.6f", tmp))
+ println(aa)
end
+ print(STDOUT, aa)
end
-function issteadystate(model::Model, oo::Output, ys::Vector{Float64})
- residuals = zeros(Float64, length(ys))
- compute_static_model_residuals!(model, oo, ys, residuals)
- return maximum(abs(residuals))<1e-6
+function display_steady_state(model::Model, oo::Output, ys::Vector{Float64})
+ aa = AutoAlign(align = Dict(1 => left, 2=> left, :default => right))
+ println("\n Deterministic Steady State:\n")
+ for i = 1:model.endo_nbr
+ if abs(ys[i])<1.0e-6
+ tmp = zero(Float64)
+ else
+ tmp = ys[i]
+ end
+ print(aa, " ", model.endo[i].name, "\t=\t", @sprintf("%.6f", tmp))
+ println(aa)
+ end
+ print(STDOUT, aa)
end
-function compute_static_model_residuals!(model::Model, oo::Output, ys::Vector{Float64}, residuals::Vector{Float64})
- model.static(ys, oo.exo_steady_state, model.params, residuals)
+function issteadystate(model::Model, oo::Output, ys::Vector{Float64})
+ r = zeros(Float64, model.endo_nbr)
+ t = zeros(Float64, model.temporaries.static[1])
+ model.static(t, r, ys, oo.exo_steady_state, model.params)
+ return norm(r)<1e-6
end
end
diff --git a/src/TestFunctions.jl b/src/TestFunctions.jl
new file mode 100644
index 0000000000000000000000000000000000000000..c8068cc46a6602f0168a57956432a5e33f3bc32e
--- /dev/null
+++ b/src/TestFunctions.jl
@@ -0,0 +1,577 @@
+module TestFunctions
+
+import Base.sign
+
+export rosenbrock, rosenbrock!, powell1, powell1!, powell2, powell2!, wood, wood!
+export helicalvalley, helicalvalley!, watson, watson!, chebyquad, chebyquad!
+export brown, brown!, discreteboundaryvalue, discreteboundaryvalue!
+export discreteintegralequation, discreteintegralequation!
+export trigonometric, trigonometric!, variablydimensioned, variablydimensioned!
+export broydentridiagonal, broydentridiagonal!
+export broydenbanded, broydenbanded!
+
+const ZERO = 0.0
+const ONE = 1.0
+const TWO = 2.0
+const THREE = 3.0
+const FOUR = 4.0
+const FIVE = 5.0
+const SIX = 6.0
+const EIGHT = 8.0
+const TEN = 10.0
+const FIFTN = 15.0
+const TWENTY = 20.0
+const HUNDRD = 100.0
+const HALF = 0.5
+const C1 = 1.0e4
+const C2 = 1.0001
+const C3 = 2.0e2
+const C4 = 20.2
+const C5 = 19.8
+const C6 = 180.0
+const C7 = 0.25
+const C8 = 0.5
+const C9 = 29.0
+const TPI = EIGHT*atan(ONE)
+
+"""
+ sign(a::Real, b::Real)
+
+Returns |a|*sign(b).
+"""
+function sign(a::Real, b::Real)
+ return abs(a)*sign(b)
+end
+
+"""
+ rosenbrock!(fval::Vector{Float64}, x::Vector{Float64})
+
+Evaluates the Rosenbrock function.
+"""
+function rosenbrock!(fval::Vector{Float64}, x::Vector{Float64})
+ fval[1] = ONE-x[1]
+ fval[2] = TEN*(x[2]-fval[1]*fval[1])
+end
+
+"""
+ rosenbrock!(fval::Matrix{Float64}, x::Vector{Float64})
+
+Evaluates the jacobian matrix of the Rosenbrock function.
+"""
+function rosenbrock!(fjac::Matrix{Float64}, x::Vector{Float64})
+ fjac[1,1] = -ONE
+ fjac[1,2] = ZERO
+ fjac[2,1] = -TWENTY*x[1]
+ fjac[2,2] = TEN
+end
+
+"""
+ rosenbrock()
+
+Provides initial guess for solving the nonlinear system of equation based on the Rosenbrock function.
+"""
+function rosenbrock()
+ return [-C1, ONE]
+end
+
+function powell1!(fval::Vector{Float64}, x::Vector{Float64})
+ fval[1] = x[1]+TEN*x[2]
+ fval[2] = sqrt(FIVE)*(x[3]-x[4])
+ tmp = x[2]-TWO*x[3]
+ fval[3] = tmp*tmp
+ tmp = x[1]-x[4]
+ fval[4] = sqrt(TEN)*tmp*tmp
+end
+
+function powell1!(fjac::Matrix{Float64}, x::Vector{Float64})
+ fill!(fjac, 0.0)
+ fjac[1,1] = ONE
+ fjac[1,2] = TEN
+ fjac[2,3] = sqrt(FIVE)
+ fjac[2,4] = -fjac[2,3]
+ fjac[3,2] = TWO*(x[2]-TWO*x[3])
+ fjac[3,3] = -TWO*fjac[3,2]
+ fjac[4,1] = TWO*sqrt(TEN)*(x[1]-x[4])
+ fjac[4,4] = -fjac[4,1]
+end
+
+function powell1()
+ return [THREE, -ONE, ZERO, ONE]
+end
+
+function powell2!(fval::Vector{Float64}, x::Vector{Float64})
+ fval[1] = C1*x[1]*x[2]-ONE
+ fval[2] = exp(-x[1])+exp(-x[2])-C2
+end
+
+function powell2!(fjac::Matrix{Float64}, x::Vector{Float64})
+ fjac[1,1] = C1*x[2]
+ fjac[1,2] = C1*x[1]
+ fjac[2,1] = -exp(-x[1])
+ fjac[2,2] = -exp(-x[2])
+end
+
+function powell2()
+ return [ZERO, ONE]
+end
+
+function wood!(fval::Vector{Float64}, x::Vector{Float64})
+ tmp1 = x[2]-x[1]*x[1]
+ tmp2 = x[4]-x[3]*x[3]
+ fval[1] = -C3*x[1]*tmp1-(ONE-x[1])
+ fval[2] = C3*tmp1+C4*(x[2]-1.0)+C5*(x[4]-ONE)
+ fval[3] = -C6*x[3]*tmp2-(ONE-x[3])
+ fval[4] = C6*tmp2+C4*(x[4]-ONE)+C5*(x[2]-ONE)
+end
+
+function wood!(fjac::Matrix{Float64}, x::Vector{Float64})
+ fill!(fjac, ZERO)
+ tmp1 = x[2]-THREE*x[1]*x[1]
+ tmp2 = x[4]-THREE*x[3]*x[3]
+ fjac[1,1] = -C3*tmp1+ONE
+ fjac[1,2] = -C3*x[1]
+ fjac[2,1] = -TWO*C3*x[1]
+ fjac[2,2] = C3+C4
+ fjac[2,4] = C5
+ fjac[3,3] = -C6*tmp2+ONE
+ fjac[3,4] = -C6*x[3]
+ fjac[4,2] = C5
+ fjac[4,3] = -TWO*C6*x[3]
+ fjac[4,4] = C6+C4
+end
+
+function wood()
+ return [-THREE, -ONE, -THREE, -ONE]
+end
+
+function helicalvalley!(fval::Vector{Float64}, x::Vector{Float64})
+ tmp1 = sign(C7, x[2])
+ if x[1]>ZERO
+ tmp1 = atan(x[2]/x[1])/TPI
+ end
+ if x[1]<ZERO
+ tmp1 = atan(x[2]/x[1])/TPI+C8
+ end
+ tmp2 = sqrt(x[1]*x[1]+x[2]*x[2])
+ fval[1] = TEN*(x[3]-TEN*tmp1)
+ fval[2] = TEN*(tmp2-ONE)
+ fval[3] = x[3]
+end
+
+function helicalvalley!(fjac::Matrix{Float64}, x::Vector{Float64})
+ tmp = x[1]*x[1]+x[2]*x[2]
+ tmp1 = TPI*tmp
+ tmp2 = sqrt(tmp)
+ fjac[1,1] = HUNDRD*x[2]/tmp1
+ fjac[1,2] = -HUNDRD*x[1]/tmp1
+ fjac[1,3] = TEN
+ fjac[2,1] = TEN*x[1]/tmp2
+ fjac[2,2] = TEN*x[2]/tmp2
+ fjac[2,3] = ZERO
+ fjac[3,1] = ZERO
+ fjac[3,2] = ZERO
+ fjac[3,3] = ONE
+end
+
+function helicalvalley()
+ return [-ONE, ZERO, ZERO]
+end
+
+function watson!(fval::Vector{Float64}, x::Vector{Float64})
+ n = length(fval)
+ fill!(fval, ZERO)
+ for i=1:29
+ ti = Float64(i)/C9
+ sum1 = ZERO
+ tmp = ONE
+ for j=2:n
+ sum1 += Float64(j-1)*tmp*x[j]
+ tmp *= ti
+ end
+ sum2 = ZERO
+ tmp = ONE
+ for j=1:n
+ sum2 += tmp*x[j]
+ tmp *= ti
+ end
+ tmp1 = sum1-sum2*sum2-ONE
+ tmp2 = TWO*ti*sum2
+ tmp = ONE/ti
+ for k=1:n
+ fval[k] += tmp*(Float64(k-1)-tmp2)*tmp1
+ tmp *= ti
+ end
+ end
+ tmp = x[2]-x[1]*x[1]-ONE
+ fval[1] += x[1]*(ONE-TWO*tmp)
+ fval[2] += tmp
+end
+
+function watson!(fjac::Matrix{Float64}, x::Vector{Float64})
+ n = length(x)
+ fill!(fjac, ZERO)
+ for i=1:29
+ ti = Float64(i)/C9
+ sum1 = ZERO
+ tmp = ONE
+ for j=2:n
+ sum1 += Float64(j-1)*tmp*x[j]
+ tmp *= ti
+ end
+ sum2 = ZERO
+ tmp = ONE
+ for j=1:n
+ sum2 += tmp*x[j]
+ tmp *= ti
+ end
+ tmp1 = TWO*(sum1-sum2*sum2-ONE)
+ tmp2 = TWO*sum2
+ tmp = ti*ti
+ tk = ONE
+ for k=1:n
+ tj = tk
+ for j=k:n
+ fjac[k,j] += tj*((Float64(k-1)/ti-tmp2)*(Float64(j-1)/ti-tmp2)-tmp1)
+ tj *= ti
+ end
+ tk *= tmp
+ end
+ end
+ fjac[1,1] = fjac[1,1]+SIX*x[1]*x[1]-TWO*x[2]+THREE
+ fjac[1,2] = fjac[1,2]-TWO*x[1]
+ fjac[2,2] = fjac[2,2]+ONE
+ for k=1:n
+ for j=k:n
+ fjac[j,k] = fjac[k,j]
+ end
+ end
+end
+
+function watson(n::Int)
+ return zeros(Float64, n)
+end
+
+function chebyquad!(fval::Vector{Float64}, x::Vector{Float64})
+ n = length(fval)
+ fill!(fval, ZERO)
+ for j=1:n
+ tmp1 = ONE
+ tmp2 = TWO*x[j]-ONE
+ tmp = TWO*tmp2
+ for i=1:n
+ fval[i] += tmp2
+ ti = tmp*tmp2-tmp1
+ tmp1 = tmp2
+ tmp2 = ti
+ end
+ end
+ tk = ONE/Float64(n)
+ iev = -1
+ for k=1:n
+ fval[k] *= tk
+ if iev>0
+ fval[k] += ONE/(Float64(k)^2-ONE)
+ end
+ iev = -iev
+ end
+end
+
+function chebyquad!(fjac::Matrix{Float64}, x::Vector{Float64})
+ n = length(x)
+ tk = ONE/Float64(n)
+ for j=1:n
+ tmp1 = ONE
+ tmp2 = TWO*x[j]-ONE
+ tmp = TWO*tmp2
+ tmp3 = ZERO
+ tmp4 = TWO
+ for k=1:n
+ fjac[k,j] = tk*tmp4
+ ti = FOUR*tmp2+tmp*tmp4-tmp3
+ tmp3 = tmp4
+ tmp4 = ti
+ ti = tmp*tmp2-tmp1
+ tmp1 = tmp2
+ Tmp2 = ti
+ end
+ end
+end
+
+function chebyquad(n::Int)
+ h = ONE/Float64(n+1)
+ x = Vector{Float64}(n)
+ for j=1:n
+ x[j] = Float64(j)*h
+ end
+ return x
+end
+
+function brown!(fval::Vector{Float64}, x::Vector{Float64})
+ n = length(fval)
+ sum = -Float64(n+1)
+ prod = ONE
+ for j=1:n
+ sum += x[j]
+ prod *= x[j]
+ end
+ for k=1:n-1
+ fval[k]=x[k]+sum
+ end
+ fval[n] = prod-ONE
+end
+
+function brown!(fjac::Matrix{Float64}, x::Vector{Float64})
+ n = length(x)
+ prod = ONE
+ for j=1:n
+ prod *= x[j]
+ for k=1:n
+ fjac[k,j] = ONE
+ end
+ fjac[j,j] = TWO
+ end
+ for j=1:n
+ tmp = x[j]
+ if abs(tmp)<eps(Float64)
+ tmp = ONE
+ prod = ONE
+ for k=1:n
+ if k!=j
+ prod *= x[k]
+ end
+ end
+ end
+ fjac[n,j] = prod/tmp
+ end
+end
+
+function brown(n::Int)
+ x = Vector{Float64}(n)
+ fill!(x, HALF)
+ return x
+end
+
+function discreteboundaryvalue!(fval::Vector{Float64}, x::Vector{Float64})
+ n = length(fval)
+ h = ONE/Float64(n+1)
+ hh = h*h
+ for k=1:n
+ tmp = (x[k]+Float64(k)*h+ONE)^3
+ tmp1 = ZERO
+ if k!=1
+ tmp1 = x[k-1]
+ end
+ tmp2 = ZERO
+ if k!=n
+ tmp2 = x[k+1]
+ end
+ fval[k] = TWO*x[k]-tmp1-tmp2+tmp*hh/TWO
+ end
+end
+
+function discreteboundaryvalue!(fjac::Matrix{Float64}, x::Vector{Float64})
+ n = length(x)
+ h = ONE/Float64(n+1)
+ hh = h*h
+ fill!(fjac, ZERO)
+ for k=1:n
+ tmp = THREE*(x[k]+Float64(k)*h+ONE)^2
+ fjac[k,k] = TWO+tmp*hh/TWO
+ if k!=1
+ fjac[k,k-1] = -ONE
+ end
+ if k!=n
+ fjac[k,k+1] = -ONE
+ end
+ end
+end
+
+function discreteboundaryvalue(n::Int)
+ h = ONE/Float64(n+1)
+ x = Vector{Float64}(n)
+ for j=1:n
+ tj = Float64(j)*h
+ x[j] = tj*(tj-ONE)
+ end
+ return x
+end
+
+function discreteintegralequation!(fval::Vector{Float64}, x::Vector{Float64})
+ n = length(fval)
+ h = ONE/Float64(n+1)
+ for k=1:n
+ tk = Float64(k)*h
+ sum1 = ZERO
+ for j=1:k
+ tj = Float64(j)*h
+ tmp = (x[j]+tj+ONE)^3
+ sum1 += tj*tmp
+ end
+ sum2 = ZERO
+ kp1 = k+1
+ if n>=kp1
+ for j=kp1:n
+ tj = Float64(j)*h
+ tmp = (x[j]+tj+ONE)^3
+ sum2 += (ONE-tj)*tmp
+ end
+ end
+ fval[k] = x[k] + h*((ONE-tk)*sum1+tk*sum2)/TWO
+ end
+end
+
+function discreteintegralequation!(fjac::Matrix{Float64}, x::Vector{Float64})
+ n = length(x)
+ h = ONE/Float64(n+1)
+ for k=1:n
+ tk = Float64(k)*h
+ for j=1:n
+ tj = Float64(j)*h
+ tmp = THREE*(x[j]+tj+ONE)^2
+ fjac[k,j] = h*min(tj*(ONE-tk),tk*(ONE-tj))*tmp/TWO
+ end
+ fjac[k,k] += ONE
+ end
+end
+
+discreteintegralequation(n::Int) = discreteboundaryvalue(n)
+
+function trigonometric!(fval::Vector{Float64}, x::Vector{Float64})
+ n = length(fval)
+ sum = ZERO
+ for j=1:n
+ fval[j] = cos(x[j])
+ sum += fval[j]
+ end
+ for k=1:n
+ fval[k] = Float64(n+k)-sin(x[k])-sum-Float64(k)*fval[k]
+ end
+end
+
+function trigonometric!(fjac::Matrix{Float64}, x::Vector{Float64})
+ n = length(x)
+ for j=1:n
+ tmp = sin(x[j])
+ for k=1:n
+ fjac[k,j] = tmp
+ end
+ fjac[j,j] = Float64(j+1)*tmp-cos(x[j])
+ end
+end
+
+function trigonometric(n::Int)
+ x = Vector{Float64}(n)
+ fill!(x, ONE/Float64(n))
+ return x
+end
+
+function variablydimensioned!(fval::Vector{Float64}, x::Vector{Float64})
+ n = length(fval)
+ sum = ZERO
+ for j=1:n
+ sum += Float64(j)*(x[j]-ONE)
+ end
+ tmp = sum*(ONE+TWO*sum*sum)
+ for k=1:n
+ fval[k] = x[k]-ONE+Float64(k)*tmp
+ end
+end
+
+function variablydimensioned!(fjac::Matrix{Float64}, x::Vector{Float64})
+ n = length(x)
+ sum = ZERO
+ for j=1:n
+ sum += Float64(j)*(x[j]-ONE)
+ end
+ tmp = ONE+SIX*sum*sum
+ for k=1:n
+ for j=1:n
+ fjac[k,j] = Float64(k*j)*tmp
+ fjac[j,k] = fjac[k,j]
+ end
+ fjac[k,k] += ONE
+ end
+end
+
+function variablydimensioned(n::Int)
+ h = ONE/Float64(n)
+ x = Vector{Float64}(n)
+ for j=1:n
+ x[j] = ONE-Float64(j)*h
+ end
+ return x
+end
+
+function broydentridiagonal!(fval::Vector{Float64}, x::Vector{Float64})
+ n = length(fval)
+ for k=1:n
+ tmp = (THREE-TWO*x[k])*x[k]
+ tmp1 = ZERO
+ if k!=1
+ tmp1 = x[k-1]
+ end
+ tmp2 = ZERO
+ if k!=n
+ tmp2 = x[k+1]
+ end
+ fval[k] = tmp-tmp1-TWO*tmp2+ONE
+ end
+end
+
+function broydentridiagonal!(fjac::Matrix{Float64}, x::Vector{Float64})
+ n = length(x)
+ fill!(fjac, ZERO)
+ for k=1:n
+ fjac[k,k] = THREE-FOUR*x[k]
+ if k!=1
+ fjac[k,k-1] = -ONE
+ end
+ if k!=n
+ fjac[k,k+1] = -TWO
+ end
+ end
+end
+
+function broydentridiagonal(n::Int)
+ x = Vector{Float64}(n)
+ fill!(x, -ONE)
+ return x
+end
+
+function broydenbanded!(fval::Vector{Float64}, x::Vector{Float64})
+ n = length(fval)
+ ml = 5
+ mu = 1
+ for k=1:n
+ k1 = max(1,k-ml)
+ k2 = min(k+mu,n)
+ tmp = ZERO
+ for j=k1:k2
+ if j!=k
+ tmp += x[j]*(ONE+x[j])
+ end
+ end
+ fval[k] = x[k]*(TWO+FIVE*x[k]*x[k])+ONE-tmp
+ end
+end
+
+function broydenbanded!(fjac::Matrix{Float64}, x::Vector{Float64})
+ n = length(x)
+ ml = 5
+ mu = 1
+ fill!(fjac,ZERO)
+ for k=1:n
+ k1 = max(1,k-ml)
+ k2 = min(k+mu,n)
+ for j=k1:k2
+ if j!=k
+ fjac[k,j] = -(ONE+TWO*x[j])
+ end
+ end
+ fjac[k,k] = TWO+FIFTN*x[k]*x[k]
+ end
+end
+
+broydenbanded(n::Int) = broydentridiagonal(n)
+
+end
diff --git a/test/dynare-solvers/test-1.jl b/test/dynare-solvers/test-1.jl
new file mode 100644
index 0000000000000000000000000000000000000000..2f691f6d5a2e1ca0dff3628fd554891bc0c5073e
--- /dev/null
+++ b/test/dynare-solvers/test-1.jl
@@ -0,0 +1,218 @@
+#=
+** This script performs tests for nonlinear solvers.
+**
+** Copyright (C) 2018 Dynare Team
+**
+** This file is part of Dynare.
+**
+** Dynare is free software: you can redistribute it and/or modify
+** it under the terms of the GNU General Public License as published by
+** the Free Software Foundation, either version 3 of the License, or
+** (at your option) any later version.
+**
+** Dynare is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+=#
+
+rootdir = @__DIR__
+origdir = pwd()
+
+if isempty(findin([abspath("$(rootdir)/../../src")], LOAD_PATH))
+ unshift!(LOAD_PATH, abspath("$(rootdir)/../../src"))
+end
+
+using Base.Test
+using TestFunctions
+
+cd("$(rootdir)")
+
+@testset "test unit-test-routines" begin
+ @test begin
+ try
+ x = rosenbrock()
+ f = Vector{Float64}(2)
+ j = Matrix{Float64}(2,2)
+ rosenbrock!(f, x)
+ rosenbrock!(j, x)
+ true
+ catch
+ false
+ end
+ end
+ @test begin
+ try
+ x = powell1()
+ n = length(x)
+ f = Vector{Float64}(n)
+ j = Matrix{Float64}(n,n)
+ powell1!(f, x)
+ powell1!(j, x)
+ true
+ catch
+ false
+ end
+ end
+ @test begin
+ try
+ x = powell2()
+ n = length(x)
+ f = Vector{Float64}(n)
+ j = Matrix{Float64}(n,n)
+ powell2!(f, x)
+ powell2!(j, x)
+ true
+ catch
+ false
+ end
+ end
+ @test begin
+ try
+ x = wood()
+ n = length(x)
+ f = Vector{Float64}(n)
+ j = Matrix{Float64}(n,n)
+ wood!(f, x)
+ wood!(j, x)
+ true
+ catch
+ false
+ end
+ end
+ @test begin
+ try
+ x = helicalvalley()
+ n = length(x)
+ f = Vector{Float64}(n)
+ j = Matrix{Float64}(n,n)
+ helicalvalley!(f, x)
+ helicalvalley!(j, x)
+ true
+ catch
+ false
+ end
+ end
+ @test begin
+ try
+ x = watson(10)
+ n = length(x)
+ f = Vector{Float64}(n)
+ j = Matrix{Float64}(n,n)
+ watson!(f, x)
+ watson!(j, x)
+ true
+ catch
+ false
+ end
+ end
+ @test begin
+ try
+ x = chebyquad(10)
+ n = length(x)
+ f = Vector{Float64}(n)
+ j = Matrix{Float64}(n,n)
+ chebyquad!(f, x)
+ chebyquad!(j, x)
+ true
+ catch
+ false
+ end
+ end
+ @test begin
+ try
+ x = brown(10)
+ n = length(x)
+ f = Vector{Float64}(n)
+ j = Matrix{Float64}(n,n)
+ brown!(f, x)
+ brown!(j, x)
+ true
+ catch
+ false
+ end
+ end
+ @test begin
+ try
+ x = discreteboundaryvalue(10)
+ n = length(x)
+ f = Vector{Float64}(n)
+ j = Matrix{Float64}(n,n)
+ discreteboundaryvalue!(f, x)
+ discreteboundaryvalue!(j, x)
+ true
+ catch
+ false
+ end
+ end
+ @test begin
+ try
+ x = discreteintegralequation(10)
+ n = length(x)
+ f = Vector{Float64}(n)
+ j = Matrix{Float64}(n,n)
+ discreteintegralequation!(f, x)
+ discreteintegralequation!(j, x)
+ true
+ catch
+ false
+ end
+ end
+ @test begin
+ try
+ x = trigonometric(10)
+ n = length(x)
+ f = Vector{Float64}(n)
+ j = Matrix{Float64}(n,n)
+ trigonometric!(f, x)
+ trigonometric!(j, x)
+ true
+ catch
+ false
+ end
+ end
+ @test begin
+ try
+ x = variablydimensioned(10)
+ n = length(x)
+ f = Vector{Float64}(n)
+ j = Matrix{Float64}(n,n)
+ variablydimensioned!(f, x)
+ variablydimensioned!(j, x)
+ true
+ catch
+ false
+ end
+ end
+ @test begin
+ try
+ x = broydentridiagonal(10)
+ n = length(x)
+ f = Vector{Float64}(n)
+ j = Matrix{Float64}(n,n)
+ broydentridiagonal!(f, x)
+ broydentridiagonal!(j, x)
+ true
+ catch
+ false
+ end
+ end
+ @test begin
+ try
+ x = broydenbanded(10)
+ n = length(x)
+ f = Vector{Float64}(n)
+ j = Matrix{Float64}(n,n)
+ broydenbanded!(f, x)
+ broydenbanded!(j, x)
+ true
+ catch
+ false
+ end
+ end
+end
+
+cd(origdir)
diff --git a/test/dynare-solvers/test-2.jl b/test/dynare-solvers/test-2.jl
new file mode 100644
index 0000000000000000000000000000000000000000..d88d95a9a5c48fb9fabfaa078177a856ed518f2f
--- /dev/null
+++ b/test/dynare-solvers/test-2.jl
@@ -0,0 +1,94 @@
+#=
+** This script performs tests for the trust region nonlinear solver.
+**
+** Copyright (C) 2018 Dynare Team
+**
+** This file is part of Dynare.
+**
+** Dynare is free software: you can redistribute it and/or modify
+** it under the terms of the GNU General Public License as published by
+** the Free Software Foundation, either version 3 of the License, or
+** (at your option) any later version.
+**
+** Dynare is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+=#
+
+rootdir = @__DIR__
+origdir = pwd()
+
+if isempty(findin([abspath("$(rootdir)/../../src")], LOAD_PATH))
+ unshift!(LOAD_PATH, abspath("$(rootdir)/../../src"))
+end
+
+using Base.Test
+using DynareSolvers
+using TestFunctions
+
+cd("$(rootdir)")
+
+@testset "test trust-region-solver" begin
+ @test begin
+ x, info = trustregion(rosenbrock!, rosenbrock!, rosenbrock(), 1.0e-6, 1.0e-6, 50)
+ info==1
+ end
+ @test begin
+ x, info = trustregion(powell1!, powell1!, powell1(), 1.0e-6, 1.0e-6, 50)
+ info==1
+ end
+ @test begin
+ x, info = trustregion(powell2!, powell2!, powell2(), 1.0e-6, 1.0e-6, 50)
+ info==1
+ end
+ @test begin
+ x, info = trustregion(wood!, wood!, wood(), 1.0e-6, 1.0e-6, 50)
+ info==1
+ end
+ @test begin
+ x, info = trustregion(helicalvalley!, helicalvalley!, helicalvalley(), 1.0e-6, 1.0e-6, 50)
+ info==1
+ end
+ @test begin
+ x, info = trustregion(watson!, watson!, watson(10), 1.0e-6, 1.0e-6, 50)
+ info==1
+ end
+ @test begin
+ x, info = trustregion(chebyquad!, chebyquad!, chebyquad(10), 1.0e-6, 1.0e-6, 50)
+ info==1
+ end
+ @test begin
+ x, info = trustregion(brown!, brown!, brown(10), 1.0e-6, 1.0e-6, 50)
+ info==1
+ end
+ @test begin
+ x, info = trustregion(discreteboundaryvalue!, discreteboundaryvalue!, discreteboundaryvalue(10), 1.0e-6, 1.0e-6, 50)
+ info==1
+ end
+ @test begin
+ x, info = trustregion(discreteintegralequation!, discreteintegralequation!, discreteintegralequation(10), 1.0e-6, 1.0e-6, 50)
+ info==1
+ end
+ @test begin
+ x, info = trustregion(trigonometric!, trigonometric!, trigonometric(10), 1.0e-6, 1.0e-6, 50)
+ info==1
+ end
+ @test begin
+ x, info = trustregion(variablydimensioned!, variablydimensioned!, variablydimensioned(10), 1.0e-6, 1.0e-6, 50)
+ info==1
+ end
+ @test begin
+ x, info = trustregion(broydentridiagonal!, broydentridiagonal!, broydentridiagonal(10), 1.0e-6, 1.0e-6, 50)
+ info==1
+ end
+ @test begin
+ x, info = trustregion(broydenbanded!, broydenbanded!, broydenbanded(10), 1.0e-6, 1.0e-6, 50)
+ info==1
+ end
+end
+
+cd(origdir)
diff --git a/test/runtests.jl b/test/runtests.jl
index a0b1146a3f2fde48348db6af44e71190b07e0eda..15295b7c8a4118f79ff1d3b4e9a339be9a6e2155 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -14,4 +14,8 @@ using Base.Test
include("$(rootdir0)/preprocessor/test-2.jl")
include("$(rootdir0)/preprocessor/test-3.jl")
include("$(rootdir0)/preprocessor/test-4.jl")
+ include("$(rootdir0)/dynare-solvers/test-1.jl")
+ include("$(rootdir0)/dynare-solvers/test-2.jl")
+ include("$(rootdir0)/steady-state/test-1.jl")
+
end
diff --git a/test/steady-state/example1.mod b/test/steady-state/example1.mod
new file mode 100644
index 0000000000000000000000000000000000000000..7d14d380d77f24c2dffd1e685622007ef84478c4
--- /dev/null
+++ b/test/steady-state/example1.mod
@@ -0,0 +1,37 @@
+// Test trust region nonlinear solver (solve_algo=4)
+var y, c, k, a, h, b;
+varexo e, u;
+
+parameters beta, rho, alpha, delta, theta, psi, tau;
+
+alpha = 0.36;
+rho = 0.95;
+tau = 0.025;
+beta = 0.99;
+delta = 0.025;
+psi = 0;
+theta = 2.95;
+
+phi = 0.1;
+
+model;
+c*theta*h^(1+psi)=(1-alpha)*y;
+k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1)))
+ *(exp(b(+1))*alpha*y(+1)+(1-delta)*k));
+y = exp(a)*(k(-1)^alpha)*(h^(1-alpha));
+k = exp(b)*(y-c)+(1-delta)*k(-1);
+a = rho*a(-1)+tau*b(-1) + e;
+b = tau*a(-1)+rho*b(-1) + u;
+end;
+
+initval;
+y = 1;
+c = 0.8;
+h = 0.3;
+k = 10;
+a = 0;
+b = 0;
+e = 0;
+u = 0;
+end;
+
diff --git a/test/steady-state/test-1.jl b/test/steady-state/test-1.jl
new file mode 100644
index 0000000000000000000000000000000000000000..8d4202f3027cb83dff5c30fd3d13ae5e2c7c410b
--- /dev/null
+++ b/test/steady-state/test-1.jl
@@ -0,0 +1,110 @@
+#=
+** This script performs tests for the steady state numerical solver.
+**
+** Copyright (C) 2018 Dynare Team
+**
+** This file is part of Dynare.
+**
+** Dynare is free software: you can redistribute it and/or modify
+** it under the terms of the GNU General Public License as published by
+** the Free Software Foundation, either version 3 of the License, or
+** (at your option) any later version.
+**
+** Dynare is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+=#
+
+rootdir = @__DIR__
+origdir = pwd()
+
+if isempty(findin([abspath("$(rootdir)/../../src")], LOAD_PATH))
+ unshift!(LOAD_PATH, abspath("$(rootdir)/../../src"))
+end
+
+using Base.Test
+using Dynare
+using SteadyState
+using DynareUnitTests
+
+cd("$(rootdir)")
+
+cp("example1.mod", "example1_1.mod"; remove_destination=true)
+
+steadystate = [1.080683, 0.803593, 11.083609, 0.000000, 0.291756, 0.000000]
+
+@testset "test steadystate-solver" begin
+ @test begin
+ try
+ @dynare "example1_1.mod"
+ true
+ catch
+ false
+ end
+ end
+ model = example1_1.model_
+ oo = example1_1.oo_
+ @test begin
+ try
+ yinit = [1.0, .8, 10.0, 0.0, .3, 0.1]
+ steady_state(model, oo, yinit, false)
+ true
+ catch
+ false
+ end
+ end
+ @test begin
+ try
+ yinit = [2.0, .8, 6.0, 0.0, .3, 0.1]
+ steady!(model, oo, yinit)
+ all((oo.steady_state-steadystate).<1.0e-6)
+ catch
+ false
+ end
+ end
+ @test begin
+ try
+ yinit = [.5, .3, 12.0, 0.0, .2, -0.1]
+ steady!(model, oo, yinit)
+ all((oo.steady_state-steadystate).<1.0e-6)
+ catch
+ false
+ end
+ end
+ @test begin
+ try
+ yinit = [.2, .3, 4.0, 0.0, .2, -0.1]
+ steady!(model, oo, yinit)
+ all((oo.steady_state-steadystate).<1.0e-6)
+ catch
+ false
+ end
+ end
+ @test begin
+ try
+ yinit = [.1, .3, 4.0, 0.0, .3, 0.1]
+ steady!(model, oo, yinit)
+ all((oo.steady_state-steadystate).<1.0e-6)
+ catch
+ false
+ end
+ end
+ @test begin
+ try
+ yinit = [20.0, .1, 30.0, 0.0, .3, 0.1]
+ steady!(model, oo, yinit)
+ all((oo.steady_state-steadystate).<1.0e-6)
+ catch
+ false
+ end
+ end
+end
+
+DynareUnitTests.clean("example1_1")
+rm("example1_1.mod")
+
+cd(origdir)