diff --git a/src/Dynare.jl b/src/Dynare.jl
index e17dfc47fbd20750c0ec9eb46f8e690ab5863299..09c36c7ad6a6f029529f43bf0f5783ba856d0bcb 100644
--- a/src/Dynare.jl
+++ b/src/Dynare.jl
@@ -1,13 +1,14 @@
module Dynare
+using LinearAlgebra
using NLsolve
-using Suppressor
+using Printf
export @var, @varexo, @parameters, @model, @modfile
const MathExpr = Union{Expr,Symbol,Number}
-type Model
+mutable struct Model
endo::Vector{Symbol}
exo::Vector{Symbol}
param::Vector{Symbol}
@@ -40,14 +41,14 @@ type Model
dynamic_mg!::Function # Dynamic model Jacobian (yy,x,p,jacobian)
end
-Model() = Model(Array{Symbol}(0), Array{Symbol}(0), Array{Symbol}(0), Array{MathExpr}(0),
+Model() = Model(Array{Symbol}(undef, 0), Array{Symbol}(undef, 0), Array{Symbol}(undef, 0), Array{MathExpr}(undef, 0),
0, 0, 0,
- Array{Int}(0), Array{Int}(0),
- Array{Int}(0), Array{Int}(0),
- Array{Int}(0), Array{Int}(0),
- Array{Int}(0),
- Array{Int}(0), Array{Int}(0),
- Array{Int}(0), Array{Int}(0),
+ Array{Int}(undef, 0), Array{Int}(undef, 0),
+ Array{Int}(undef, 0), Array{Int}(undef, 0),
+ Array{Int}(undef, 0), Array{Int}(undef, 0),
+ Array{Int}(undef, 0),
+ Array{Int}(undef, 0), Array{Int}(undef, 0),
+ Array{Int}(undef, 0), Array{Int}(undef, 0),
0, 0, 0, 0, 0, 0, 0,
(y::Vector{Float64}, x::Vector{Float64}, p::Vector{Float64}, fvec::Vector{Float64}) -> nothing,
(y::Vector{Float64}, x::Vector{Float64}, p::Vector{Float64}, fjac::Matrix{Float64}) -> nothing,
@@ -71,10 +72,13 @@ end
remove_line_blocks(e::Any) = e
function remove_line_blocks(e::Expr)
- if e.head == :block
- assert(sum(map(y->(typeof(y) == Symbol || y.head != :line), e.args)) == 1)
+ if typeof(e) == LineNumberNode
+ return
+ elseif e.head == :block
+# assert(sum(map(y->(typeof(y) == Symbol || y.head != :line), e.args)) == 1)
for x in e.args
- if typeof(x) == Symbol || x.head != :line
+# if typeof(x) == Symbol || x.head != :line
+ if typeof(x) != LineNumberNode
return x
end
end
@@ -99,7 +103,7 @@ scan_expr(e::Number, endo::Vector{Symbol}, lead_lags::Dict{Symbol,Tuple{Int,Int}
scan_expr(e::Symbol, endo::Vector{Symbol}, lead_lags::Dict{Symbol,Tuple{Int,Int}}) = nothing
function scan_expr(e::Expr, endo::Vector{Symbol}, lag_lead::Dict{Symbol,Tuple{Int,Int}})
- if e.head == :call
+ if typeof(e) != LineNumberNode && e.head == :call
if !isempty(filter(x->x==e.args[1], endo))
v = e.args[1]
shift = e.args[2]
@@ -175,8 +179,8 @@ import Calculus.differentiate
include("symbolic.jl")
function compute_static_mf(m::Model)
- reqs = Array{MathExpr}(m.n_endo)
- jacob = Array{MathExpr}(m.n_endo, m.n_endo)
+ reqs = Array{MathExpr}(undef, m.n_endo)
+ jacob = Array{MathExpr}(undef, m.n_endo, m.n_endo)
# Compute expressions for static residuals and jacobian
for i = 1:m.n_endo
@@ -205,8 +209,8 @@ function compute_static_mf(m::Model)
end
# Construct static model functions
- f_assigns = Array{Expr}(0)
- g_assigns = Array{Expr}(0)
+ f_assigns = Array{Expr}(undef, 0)
+ g_assigns = Array{Expr}(undef, 0)
for i = 1:m.n_endo
push!(f_assigns, :(fvec[$i] = $(reqs[i])))
for j = 1:m.n_endo
@@ -215,13 +219,13 @@ function compute_static_mf(m::Model)
end
@eval begin
- @suppress function static_mf!(y::Vector{Float64}, x::Vector{Float64}, p::Vector{Float64}, fvec::Vector{Float64})
+ function static_mf!(y::Vector{Float64}, x::Vector{Float64}, p::Vector{Float64}, fvec::Vector{Float64})
$(f_assigns...)
end
end
@eval begin
- @suppress function static_mg!(y::Vector{Float64}, x::Vector{Float64}, p::Vector{Float64}, fjac::Matrix{Float64})
+ function static_mg!(y::Vector{Float64}, x::Vector{Float64}, p::Vector{Float64}, fjac::Matrix{Float64})
$(g_assigns...)
end
end
@@ -233,17 +237,17 @@ end
function compute_dynamic_mf(m::Model)
ndyn = m.n_back_mixed + m.n_endo + m.n_fwrd_mixed
njcols = ndyn + m.n_exo
- reqs = Array{MathExpr}(m.n_endo)
- jacob = Array{MathExpr}(m.n_endo, njcols)
+ reqs = Array{MathExpr}(undef, m.n_endo)
+ jacob = Array{MathExpr}(undef, m.n_endo, njcols)
# Create gensyms for lagged and leaded variables
- lag_gensyms = Array{Symbol}(m.n_back_mixed)
+ lag_gensyms = Array{Symbol}(undef, m.n_back_mixed)
lag_dict = Dict{Symbol,Symbol}()
for i = 1:m.n_back_mixed
lag_gensyms[i] = gensym()
lag_dict[m.endo[m.zeta_back_mixed[i]]] = lag_gensyms[i]
end
- lead_gensyms = Array{Symbol}(m.n_fwrd_mixed)
+ lead_gensyms = Array{Symbol}(undef, m.n_fwrd_mixed)
lead_dict = Dict{Symbol,Symbol}()
for i = 1:m.n_fwrd_mixed
lead_gensyms[i] = gensym()
@@ -282,8 +286,8 @@ function compute_dynamic_mf(m::Model)
end
# Construct dynamic model functions
- f_assigns = Array{Expr}(0)
- g_assigns = Array{Expr}(0)
+ f_assigns = Array{Expr}(undef, 0)
+ g_assigns = Array{Expr}(undef, 0)
for i = 1:m.n_endo
push!(f_assigns, :(fvec[$i] = $(reqs[i])))
for j = 1:njcols
@@ -292,13 +296,13 @@ function compute_dynamic_mf(m::Model)
end
@eval begin
- @suppress function dynamic_mf!(yy::Vector{Float64}, x::Vector{Float64}, p::Vector{Float64}, fvec::Vector{Float64})
+ function dynamic_mf!(yy::Vector{Float64}, x::Vector{Float64}, p::Vector{Float64}, fvec::Vector{Float64})
$(f_assigns...)
end
end
@eval begin
- @suppress function dynamic_mg!(yy::Vector{Float64}, x::Vector{Float64}, p::Vector{Float64}, fjac::Matrix{Float64})
+ function dynamic_mg!(yy::Vector{Float64}, x::Vector{Float64}, p::Vector{Float64}, fjac::Matrix{Float64})
$(g_assigns...)
end
end
@@ -321,7 +325,7 @@ function compute_model_info(m::Model)
end
function calib2vec(m::Model, calib::Dict{Symbol, Float64})
- p = Array{Float64}(m.n_param)
+ p = Array{Float64}(undef, m.n_param)
for i in 1:m.n_param
p[i] = calib[m.param[i]]
end
@@ -329,7 +333,7 @@ function calib2vec(m::Model, calib::Dict{Symbol, Float64})
end
function initval2vec(m::Model, initval::Dict{Symbol, Float64})
- y = Array{Float64}(m.n_endo)
+ y = Array{Float64}(undef, m.n_endo)
for i in 1:m.n_endo
y[i] = initval[m.endo[i]]
end
@@ -337,7 +341,7 @@ function initval2vec(m::Model, initval::Dict{Symbol, Float64})
end
function exoval2vec(m::Model, exoval::Dict{Symbol, Float64})
- x = Array{Float64}(m.n_exo)
+ x = Array{Float64}(undef, m.n_exo)
for i in 1:m.n_exo
x[i] = exoval[m.exo[i]]
end
@@ -348,8 +352,9 @@ function steady_state(m::Model, calib::Dict{Symbol, Float64}, initval::Dict{Symb
p = calib2vec(m, calib)
iv = initval2vec(m, initval)
ev = exoval2vec(m, exoval)
- mf!(y::Vector{Float64}, fvec::Vector{Float64}) = m.static_mf!(y, ev, p, fvec)
- mg!(y::Vector{Float64}, fjac::Matrix{Float64}) = m.static_mg!(y, ev, p, fjac)
+ mf!(fvec::Vector{Float64}, y::Vector{Float64}) = m.static_mf!(y, ev, p, fvec)
+ mg!(fjac::Matrix{Float64}, y::Vector{Float64}, ) = m.static_mg!(y, ev, p, fjac)
+
r = nlsolve(mf!, mg!, iv)
@assert converged(r)
return(r.zero)
diff --git a/src/decision_rules.jl b/src/decision_rules.jl
index f080a265f7c94c5932e9af35018d10c0cc19e665..68f436d11c8827d78bdaf72c95e278a32cdec708 100644
--- a/src/decision_rules.jl
+++ b/src/decision_rules.jl
@@ -4,28 +4,28 @@ function decision_rules(m::Model, calib::Dict{Symbol, Float64}, steady_state::Ve
yy = [ steady_state[m.zeta_back_mixed]; steady_state; steady_state[m.zeta_fwrd_mixed] ]
x = zeros(m.n_exo)
p = calib2vec(m, calib)
- r = Array{Float64}(m.n_endo)
+ r = Array{Float64}(undef, m.n_endo)
m.dynamic_mf!(yy, x, p, r)
- jacob = Array{Float64}(m.n_endo, n_dynvars + m.n_exo)
+ jacob = Array{Float64}(undef, m.n_endo, n_dynvars + m.n_exo)
m.dynamic_mg!(yy, x, p, jacob)
A = jacob[:, 1:n_dynvars]
# Deal with static variables
if m.n_static > 0
- S = jacob[:, m.n_back_mixed + m.zeta_static]
- F = qrfact(S)
- @assert rank(F[:R]) == m.n_static
- A = F[:Q]'*A
+ S = jacob[:, m.n_back_mixed .+ m.zeta_static]
+ F = qr(S)
+ @assert rank(F.R) == m.n_static
+ A = F.Q'*A
end
- D = [ A[m.n_static+1:end, m.n_back_mixed + m.zeta_back_mixed] A[m.n_static+1:end, m.n_back_mixed+m.n_endo+1:end];
- eye(m.n_back_mixed, m.n_back_mixed)[m.beta_back, :] zeros(m.n_mixed, m.n_fwrd_mixed) ]
+ D = [ A[m.n_static+1:end, m.n_back_mixed .+ m.zeta_back_mixed] A[m.n_static+1:end, m.n_back_mixed+m.n_endo+1:end];
+ Matrix(I, m.n_back_mixed, m.n_back_mixed)[m.beta_back, :] zeros(m.n_mixed, m.n_fwrd_mixed) ]
- Atilde0plus = A[m.n_static+1:end, m.n_back_mixed + m.zeta_fwrd_mixed]
+ Atilde0plus = A[m.n_static + 1:end, m.n_back_mixed .+ m.zeta_fwrd_mixed]
Atilde0plus[:, m.beta_fwrd] = zeros(m.n_dynamic, m.n_mixed)
E = [ -A[m.n_static+1:end, 1:m.n_back_mixed] -Atilde0plus;
- zeros(m.n_mixed, m.n_back_mixed) eye(m.n_fwrd_mixed, m.n_fwrd_mixed)[m.beta_back, :] ]
+ zeros(m.n_mixed, m.n_back_mixed) Matrix(I, m.n_fwrd_mixed, m.n_fwrd_mixed)[m.beta_back, :] ]
(Z, sdim, eigs) = qz!(E, D)
@@ -44,22 +44,22 @@ function decision_rules(m::Model, calib::Dict{Symbol, Float64}, steady_state::Ve
S11 = E[1:m.n_back_mixed, 1:m.n_back_mixed]
gyminus = Z11p*(T11\S11)/Z11p
- gy = Array{Float64}(m.n_endo, m.n_back_mixed)
+ gy = Array{Float64}(undef, m.n_endo, m.n_back_mixed)
gy[m.zeta_back_mixed, :] = gyminus
gy[m.zeta_fwrd_mixed, :] = gyplus
-
+
if m.n_static > 0
- Abreve0s = A[1:m.n_static, m.n_back_mixed + m.zeta_static]
- Abreve0d = A[1:m.n_static, m.n_back_mixed + m.zeta_dynamic]
+ Abreve0s = A[1:m.n_static, m.n_back_mixed .+ m.zeta_static]
+ Abreve0d = A[1:m.n_static, m.n_back_mixed .+ m.zeta_dynamic]
Abreveminus = A[1:m.n_static, 1:m.n_back_mixed]
- Abreveplus = A[1:m.n_static, m.n_back_mixed+m.n_endo+1:end]
- gy[m.zeta_static, :] = -Abreve0s \ (Abreveplus*gyplus*gyminus + Abreve0d*gy[m.zeta_dynamic, :] + Abreveminus)
+ Abreveplus = A[1:m.n_static, m.n_back_mixed .+ m.n_endo+1:end]
+ gy[m.zeta_static, :] = -Abreve0s \ (Abreveplus*gyplus*gyminus .+ Abreve0d*gy[m.zeta_dynamic, :] + Abreveminus)
end
# Compute DR wrt shocks
- B = jacob[:, m.n_back_mixed + (1:m.n_endo)]
- B[:, m.zeta_back_mixed] += jacob[:, m.n_back_mixed+m.n_endo+(1:m.n_fwrd_mixed)]*gyplus
- gu = -B\jacob[:,m.n_back_mixed+m.n_endo+m.n_fwrd_mixed+1:end]
+ B = jacob[:, m.n_back_mixed .+ (1:m.n_endo)]
+ B[:, m.zeta_back_mixed] += jacob[:, m.n_back_mixed + m.n_endo .+ (1:m.n_fwrd_mixed)]*gyplus
+ gu = -B\jacob[:,m.n_back_mixed+m.n_endo+m.n_fwrd_mixed + 1:end]
(gy, gu, eigs)
end
diff --git a/src/perfect_foresight_simul.jl b/src/perfect_foresight_simul.jl
index 632af57b96fc0c5977b20381dee94b8abaca00b0..14f396f2ee3912306c83bf02b0467ce909a42538 100644
--- a/src/perfect_foresight_simul.jl
+++ b/src/perfect_foresight_simul.jl
@@ -1,3 +1,5 @@
+using SparseArrays
+
# Computes the perfect foresight simulation of a model for periods 1...T
# - endopath is a matrix with time in columns (T+2) and with endogenous variables in rows
# first and last column are respectively initial and terminal conditions
@@ -20,12 +22,12 @@ function perfect_foresight_simul!(m::Model, endopath::Matrix{Float64}, exopath::
Y = vec(endopath)
- d1 = Array{Float64}(m.n_endo)
- jacob = Array{Float64}(m.n_endo, m.n_back_mixed + m.n_endo + m.n_fwrd_mixed + m.n_exo)
+ d1 = Array{Float64}(undef, m.n_endo)
+ jacob = Array{Float64}(undef, m.n_endo, m.n_back_mixed + m.n_endo + m.n_fwrd_mixed + m.n_exo)
for iter = 1:maxit
- i_cols = [m.zeta_back_mixed; m.n_endo+(1:m.n_endo); 2*m.n_endo+m.zeta_fwrd_mixed]
+ i_cols = [m.zeta_back_mixed; m.n_endo .+ (1:m.n_endo); 2*m.n_endo .+ m.zeta_fwrd_mixed]
i_rows = 1:m.n_endo
for it = 2:T+1
@@ -35,19 +37,19 @@ function perfect_foresight_simul!(m::Model, endopath::Matrix{Float64}, exopath::
if it == T+1 && it == 2
idx = m.n_back_mixed+(1:m.n_endo)
elseif it == T+1
- idx = 1:(m.n_back_mixed+m.n_endo)
+ idx = 1:(m.n_back_mixed .+ m.n_endo)
elseif it == 2
- idx = m.n_back_mixed+1:(m.n_endo+m.n_fwrd_mixed)
+ idx = m.n_back_mixed+1:(m.n_endo .+ m.n_fwrd_mixed)
else
- idx = 1:(m.n_back_mixed+m.n_endo+m.n_fwrd_mixed)
+ idx = 1:(m.n_back_mixed+m.n_endo .+ m.n_fwrd_mixed)
end
- A[i_rows,i_cols[idx]-m.n_endo] = jacob[:, idx]
+ A[i_rows,i_cols[idx] .- m.n_endo] = jacob[:, idx]
res[i_rows] = d1
- i_rows += m.n_endo
- i_cols += m.n_endo
+ i_rows = i_rows .+ m.n_endo
+ i_cols = i_cols .+ m.n_endo
end
err = maximum(abs, res)
@@ -57,7 +59,7 @@ function perfect_foresight_simul!(m::Model, endopath::Matrix{Float64}, exopath::
return
end
- Y[m.n_endo + (1:T*m.n_endo)] += -A\res
+ Y[m.n_endo .+ (1:T*m.n_endo)] += -A\res
end
error("Convergence failed")
diff --git a/src/qz.jl b/src/qz.jl
index 87a9032af625501ecf12bcd26337b9aadd84bfef..66332466fbfe710b1a315b41cd394a8e327721b0 100644
--- a/src/qz.jl
+++ b/src/qz.jl
@@ -1,54 +1,54 @@
-import Base.LinAlg.BlasInt
-import Base.LinAlg.BLAS.@blasfunc
+import LinearAlgebra: BlasInt
+import LinearAlgebra.BLAS.@blasfunc
const liblapack = Base.liblapack_name
# Equivalent of mjdgges
# S and T are modified in place
+function selctg(alphar::Ptr{Float64}, alphai::Ptr{Float64}, beta::Ptr{Float64})
+ criterium::Float64 = 1.0 + 1e-6
+ ai = unsafe_load(alphai, 1)
+ ar = unsafe_load(alphar, 1)
+ b = unsafe_load(beta, 1)
+
+ convert(BlasInt, ar^2 + ai^2 < criterium * b^2)
+end
+
@eval begin
function qz!(S::Matrix{Float64}, T::Matrix{Float64})
- function selctg(alphar::Ptr{Float64}, alphai::Ptr{Float64}, beta::Ptr{Float64})
- const criterium::Float64 = 1.0 + 1e-6
- ai = unsafe_load(alphai, 1)
- ar = unsafe_load(alphar, 1)
- b = unsafe_load(beta, 1)
-
- convert(BlasInt, ar^2 + ai^2 < criterium * b^2)
- end
-
- selctg2 = cfunction(selctg, BlasInt, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}))
+ selctg2 = @cfunction(selctg, BlasInt, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}))
(n,) = size(S)
@assert size(S) == (n,n)
@assert size(T) == (n,n)
@assert stride(S, 1) == 1
@assert stride(T, 1) == 1
- alphar = Array{Float64}(n)
- alphai = Array{Float64}(n)
- beta = Array{Float64}(n)
+ alphar = Array{Float64}(undef, n)
+ alphai = Array{Float64}(undef, n)
+ beta = Array{Float64}(undef, n)
lwork = 16*n+16 # Same heuristic choice than mjdgges
- work = Array{Float64}(lwork)
- vsl = Array{Float64}(n, n)
- vsr = Array{Float64}(n, n)
- bwork = Array{BlasInt}(n)
- sdim = Array{BlasInt}(1)
- info = Array{BlasInt}(1)
+ work = Array{Float64}(undef, lwork)
+ vsl = Array{Float64}(undef, n, n)
+ vsr = Array{Float64}(undef, n, n)
+ bwork = Array{BlasInt}(undef, n)
+ sdim = Array{BlasInt}(undef, 1)
+ info = Array{BlasInt}(undef, 1)
- ccall((@blasfunc(dgges_), liblapack), Void, (Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{Void},
- Ptr{BlasInt}, Ptr{Float64}, Ptr{BlasInt},
- Ptr{Float64}, Ptr{BlasInt}, Ptr{BlasInt},
+ ccall((@blasfunc(dgges_), liblapack), Nothing, (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ptr{Nothing},
+ Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
+ Ptr{Float64}, Ref{BlasInt}, Ref{BlasInt},
Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
- Ptr{Float64}, Ptr{BlasInt},
- Ptr{Float64}, Ptr{BlasInt},
- Ptr{Float64}, Ptr{BlasInt},
- Ptr{BlasInt}, Ptr{BlasInt}),
- &'N', &'V', &'S', selctg2, &n, S, &stride(S, 2), T, &stride(T, 2),
- sdim, alphar, alphai, beta, vsl, &stride(vsl, 2),
- vsr, &stride(vsr, 2), work, &lwork, bwork, info)
+ Ptr{Float64}, Ref{BlasInt},
+ Ptr{Float64}, Ref{BlasInt},
+ Ptr{Float64}, Ref{BlasInt},
+ Ref{BlasInt}, Ptr{BlasInt}),
+ 'N', 'V', 'S', selctg2, n, S, stride(S, 2), T, stride(T, 2),
+ sdim, alphar, alphai, beta, vsl, stride(vsl, 2),
+ vsr, stride(vsr, 2), work, lwork, bwork, info)
@assert info[1] == 0
- eigs = Array{Complex128}(n)
+ eigs = Array{ComplexF64}(undef, n)
for i = 1:n
eigs[i] = alphar[i]/beta[i] + (alphai[i] == 0 && beta[i] == 0 ? 0.0 : alphai[i]/beta[i])*im
end
diff --git a/test/ramst.jl b/test/ramst.jl
index c7354b06d1deb28c31e888576d93f73ab467389e..52833166cf506d0342832e6017d0c7d7c6a65c2f 100644
--- a/test/ramst.jl
+++ b/test/ramst.jl
@@ -48,7 +48,7 @@ print_steady_state(m, s)
# condition, i.e. it goes from t=0 to t=201)
T = 200
-endopath = repmat(s, 1, T+2)
+endopath = repeat(s, 1, T+2)
exopath = ones(m.n_exo, T)
exopath[1, 1] = 1.2
diff --git a/test/runtests.jl b/test/runtests.jl
index 74f6ccb947c1fc2ecd0ac9eae9c0065a939b1902..10abdda719fc75a69e910a8baebf59953ef80c68 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,5 @@
-if isempty(findin([abspath("../src")], LOAD_PATH))
- unshift!(LOAD_PATH, abspath("../src"))
+if isempty(findall((in)(LOAD_PATH), abspath("../src") ))
+ push!(LOAD_PATH, abspath("../src"))
end
using Dynare