Commit c036da91 authored by MichelJuillard's avatar MichelJuillard
Browse files

adding histval and deterministic_trends

parent 5805e6b0
This diff is collapsed.
......@@ -4,12 +4,21 @@ authors = ["michel "]
version = "0.1.0"
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
FastLapackInterface = "29a986be-02c6-4525-aec4-84b980013641"
GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
Gadfly = "c91e804a-d5a3-530f-b6f0-dfbca275c004"
JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
KalmanFilterTools = "ec016732-9331-11e9-305f-d53cb440b85a"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
LinearRationalExpectations = "45f42fbc-210c-4ecc-9452-59ec793b9bfd"
Periods = "0c1d1a06-9168-427d-90e3-b4d73349562f"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PolynomialMatrixEquations = "4f9d485d-518f-41ed-81c8-372cd804c93b"
PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
TimeDataFrames = "ff13af16-2f41-468a-a2e1-aff5c43aabf7"
[compat]
PrettyTables = "0.8.4"
using LinearAlgebra
function simul_first_order!(results, initial_values, x, c, A, B, periods)
y_1 = initial_values
for t = 1:periods
r_1 = view(results, 1, :)
r_1 .= initial_values .- c
for t = 2:periods + 1
r = view(results, t, :)
e = view(x, t, :)
mul!(r, B, e)
mul!(r, A, y_1, 1.0, 1.0)
r .+= c
mul!(r, A, r_1, 1.0, 1.0)
r_1 .+= c
r_1 = r
end
r_1 .+= c
end
......@@ -25,6 +25,7 @@ struct Model
i_fwrd_b
i_fwrd_ns
i_both
i_non_states
p_static
p_bkwrd
p_bkwrd_b
......@@ -59,12 +60,37 @@ struct Model
i_lagged_exogenous
serially_correlated_exogenous
Sigma_e
maximum_endo_lag
maximum_endo_lead
maximum_exo_lag
maximum_exo_lead
maximum_exo_det_lag
maximum_exo_det_lead
maximum_lag
maximum_lead
orig_maximum_endo_lag
orig_maximum_endo_lead
orig_maximum_exo_lag
orig_maximum_exo_lead
orig_maximum_exo_det_lag
orig_maximum_exo_det_lead
orig_maximum_lag
orig_maximum_lead
dynamic!
static!
steady_state!
end
function Model(modfilename, endo_nbr, lead_lag_incidence, exogenous_nbr, lagged_exogenous_nbr, exogenous_deterministic_nbr, parameter_nbr)
function Model(modfilename, endo_nbr, lead_lag_incidence,
exogenous_nbr, lagged_exogenous_nbr, exogenous_deterministic_nbr,
parameter_nbr, maximum_endo_lag, maximum_endo_lead,
maximum_exo_lag, maximum_exo_lead, maximum_exo_det_lag,
maximum_exo_det_lead, maximum_lag, maximum_lead,
orig_maximum_endo_lag, orig_maximum_endo_lead,
orig_maximum_exo_lag, orig_maximum_exo_lead,
orig_maximum_exo_det_lag, orig_maximum_exo_det_lead,
orig_maximum_lag, orig_maximum_lead )
i_static = findall((lead_lag_incidence[1,:] .== 0) .& (lead_lag_incidence[3,:] .== 0))
p_static = lead_lag_incidence[2,i_static]
i_dyn = findall((lead_lag_incidence[1,:] .> 0) .| (lead_lag_incidence[3,:] .> 0))
......@@ -81,7 +107,8 @@ function Model(modfilename, endo_nbr, lead_lag_incidence, exogenous_nbr, lagged_
p_fwrd = lead_lag_incidence[3,i_fwrd]
p_fwrd_b = lead_lag_incidence[3,i_dyn[i_fwrd_ns]]
n_fwrd = length(i_fwrd)
i_both = findall((lead_lag_incidence[1,:] .> 0) .& (lead_lag_incidence[3,:] .> 0))
i_both = findall((lead_lag_incidence[1,:] .> 0) .& (lead_lag_incidence[3,:] .> 0))
i_non_states = union(i_fwrd, i_static)
p_both_b = lead_lag_incidence[1,i_both]
p_both_f = lead_lag_incidence[3,i_both]
n_both = length(i_both)
......@@ -125,19 +152,26 @@ function Model(modfilename, endo_nbr, lead_lag_incidence, exogenous_nbr, lagged_
else
steady_state! = nothing
end
Model(endo_nbr, exogenous_nbr, lagged_exogenous_nbr, exogenous_deterministic_nbr,
parameter_nbr, lead_lag_incidence, n_static, n_fwrd, n_bkwrd, n_both,
Model(endo_nbr, exogenous_nbr, lagged_exogenous_nbr,
exogenous_deterministic_nbr, parameter_nbr,
lead_lag_incidence, n_static, n_fwrd, n_bkwrd, n_both,
n_states, DErows1, DErows2, n_dyn, i_static, i_dyn, i_bkwrd,
i_bkwrd_b, i_bkwrd_ns, i_fwrd, i_fwrd_b, i_fwrd_ns, i_both,
p_static, p_bkwrd, p_bkwrd_b, p_fwrd, p_fwrd_b, p_both_b,
p_both_f, i_current, p_current, n_current, i_current_ns,
p_current_ns, n_current_ns, icolsD, jcolsD, icolsE, jcolsE,
colsUD, colsUE, i_cur_fwrd, n_cur_fwrd, p_cur_fwrd,
i_cur_bkwrd, n_cur_bkwrd, p_cur_bkwrd, i_cur_both,
n_cur_both, p_cur_both, gx_rows, hx_rows,
i_current_exogenous, i_lagged_exogenous,
serially_correlated_exogenous, Sigma_e, dynamic!, static!,
steady_state!)
i_non_states, p_static, p_bkwrd, p_bkwrd_b, p_fwrd,
p_fwrd_b, p_both_b, p_both_f, i_current, p_current,
n_current, i_current_ns, p_current_ns, n_current_ns, icolsD,
jcolsD, icolsE, jcolsE, colsUD, colsUE, i_cur_fwrd,
n_cur_fwrd, p_cur_fwrd, i_cur_bkwrd, n_cur_bkwrd,
p_cur_bkwrd, i_cur_both, n_cur_both, p_cur_both, gx_rows,
hx_rows, i_current_exogenous, i_lagged_exogenous,
serially_correlated_exogenous, Sigma_e, maximum_endo_lag,
maximum_endo_lead, maximum_exo_lag, maximum_exo_lead,
maximum_exo_det_lag, maximum_exo_det_lead, maximum_lag,
maximum_lead, orig_maximum_endo_lag, orig_maximum_endo_lead,
orig_maximum_exo_lag, orig_maximum_exo_lead,
orig_maximum_exo_det_lag, orig_maximum_exo_det_lead,
orig_maximum_lag, orig_maximum_lead, dynamic!, static!,
steady_state!)
end
function inverse_order_of_dynare_decision_rule(m::Model)
......
using CSV
using DataFrames
using GR
using JSON
using KalmanFilterTools
using LinearRationalExpectations
using Periods
using Plots
using TimeDataFrames
struct Simulation
name::String
statement::String
options::Dict{String, Any}
data::TimeDataFrame
end
struct ModelResults
endogenous_steady_state::Vector{Float64}
endogenous_trend::Matrix{Float64}
endogenous_variance::Matrix{Float64}
exogenous_steady_state::Vector{Float64}
exogenous_deterministic_steady_state::Vector{Float64}
linearrationalexpectations::LinearRationalExpectationsResults
simulations::Vector{Simulation}
smoother::Dict{String, Any}
end
struct Results
......@@ -19,6 +37,8 @@ mutable struct Work
dynamic_variables::Vector{Float64}
jacobian::Matrix{Float64}
qr_jacobian::Matrix{Float64}
model_has_trend::Bool
histval::Matrix{Float64}
end
mutable struct Context
......@@ -38,8 +58,39 @@ struct ModelInfo
nsfwrd
nspred
ndynamic
maximum_endo_lag
maximum_endo_lead
maximum_exo_lag
maximum_exo_lead
maximum_exo_det_lag
maximum_exo_det_lead
maximum_lag
maximum_lead
orig_maximum_endo_lag
orig_maximum_endo_lead
orig_maximum_exo_lag
orig_maximum_exo_lead
orig_maximum_exo_det_lag
orig_maximum_exo_det_lead
orig_maximum_lag
orig_maximum_lead
end
context = Context(Dict{String, DynareSymbol}(),
Vector{Model}(undef, 0),
Dict(),
Results(Vector{ModelResults}(undef, 0)),
Work(Vector{Float64}(undef, 0),
Vector{Float64}(undef, 0),
Vector{Float64}(undef, 0),
Vector{Float64}(undef, 0),
Matrix{Float64}(undef, 0, 0),
Matrix{Float64}(undef, 0, 0),
false,
Matrix{Float64}(undef, 0, 0)
)
)
function parser(modfilename)
modelstring = open(f -> read(f, String), modfilename*"/model/json/modfile.json")
modeljson = JSON.parse(modelstring)
......@@ -58,15 +109,53 @@ function parser(modfilename)
exo_nbr,
0,
exo_det_nbr,
param_nbr)
param_nbr,
model_info.maximum_endo_lag,
model_info.maximum_endo_lead,
model_info.maximum_exo_lag,
model_info.maximum_exo_lead,
model_info.maximum_exo_det_lag,
model_info.maximum_exo_det_lead,
model_info.maximum_lag,
model_info.maximum_lead,
model_info.orig_maximum_endo_lag,
model_info.orig_maximum_endo_lead,
model_info.orig_maximum_exo_lag,
model_info.orig_maximum_exo_lead,
model_info.orig_maximum_exo_det_lag,
model_info.orig_maximum_exo_det_lead,
model_info.orig_maximum_lag,
model_info.orig_maximum_lead
)
if "varobs" in keys(modeljson)
varobs = modeljson["varobs"]
varobs_ids = modeljson["varobs_ids"]
else
varobs = []
varobs_ids = []
end
if "varexobs" in keys(modeljson)
varexobs = modeljson["varexobs"]
varexobs_ids = modeljson["varexobs_ids"]
else
varexobs = []
varexobs_ids = []
end
order = 1
modelresults = ModelResults(Vector{Float64}(undef, endo_nbr),
Matrix{Float64}(undef, endo_nbr, 2),
Matrix{Float64}(undef, endo_nbr, endo_nbr),
Vector{Float64}(undef, exo_nbr),
Vector{Float64}(undef, exo_det_nbr),
LinearRationalExpectationsResults(order,
endo_nbr,
exo_nbr,
model.n_states))
model.n_states),
Vector{Simulation}(undef, 0),
Dict{String, Any}())
ncol = model.n_bkwrd + model.n_current + model.n_fwrd + 2*model.n_both
ncol1 = ncol + model.exogenous_nbr
work = Work(Vector{Float64}(undef, model.parameter_nbr),
......@@ -74,28 +163,39 @@ function parser(modfilename)
Vector{Float64}(undef, sum(model.dynamic!.tmp_nbr[1:2])),
Vector{Float64}(undef, ncol),
Matrix{Float64}(undef, model.endogenous_nbr, ncol1),
Matrix{Float64}(undef, model.endogenous_nbr, ncol1))
Matrix{Float64}(undef, model.endogenous_nbr, ncol1),
false,
Matrix{Float64}(undef, 0, 0))
results = Results([modelresults])
context = Context(symboltable, [model], Dict(), results, work)
global context = Context(symboltable, [model], Dict(), results, work)
for field in modeljson["statements"]
if field["statementName"] == "param_init"
initialize_parameter!(work.params, field, symboltable)
elseif field["statementName"] == "native"
native_statement(field)
elseif field["statementName"] == "initval"
initval(field)
elseif field["statementName"] == "shocks"
shocks!(model.Sigma_e, field, symboltable)
elseif field["statementName"] == "verbatim"
verbatim(field)
if field["statementName"] == "calib_smoother"
calib_smoother!(context, field, varobs, varobs_ids)
elseif field["statementName"] == "check"
check(field)
elseif field["statementName"] == "stoch_simul"
stoch_simul!(context, field)
elseif field["statementName"] == "deterministic_trends"
deterministic_trends!(context, field)
elseif field["statementName"] == "histval"
histval!(context, field, symboltable)
elseif field["statementName"] == "initval"
initval!(field)
elseif field["statementName"] == "native"
@show field["string"]
expr = Meta.parse(field["string"])
eval(expr)
native_statement(field)
elseif field["statementName"] == "param_init"
initialize_parameter!(work.params, field, symboltable)
elseif field["statementName"] == "perfect_foresight_setup"
perfect_foresight_setup!(options, field)
elseif field["statementName"] == "perfect_foresight_solver"
perfect_foresight_solver!(context, field)
elseif field["statementName"] == "shocks"
shocks!(model.Sigma_e, field, symboltable)
elseif field["statementName"] == "stoch_simul"
stoch_simul!(context, field)
elseif field["statementName"] == "verbatim"
verbatim(field)
else
error("Unrecognized statement $(field["statementName"])")
end
......@@ -119,15 +219,35 @@ function set_symbol_table!(table::Dict{String, DynareSymbol},
return count
end
get_model_info(field) = ModelInfo(hcat(field["lead_lag_incidence"]...),
field["nstatic"],
field["nfwrd"],
field["npred"],
field["nboth"],
field["nsfwrd"],
field["nspred"],
field["ndynamic"])
get_model_info(field) =
ModelInfo(hcat(field["lead_lag_incidence"]...),
field["nstatic"],
field["nfwrd"],
field["npred"],
field["nboth"],
field["nsfwrd"],
field["nspred"],
field["ndynamic"],
field["maximum_endo_lag"],
field["maximum_endo_lead"],
field["maximum_exo_lag"],
field["maximum_exo_lead"],
field["maximum_exo_det_lag"],
field["maximum_exo_det_lead"],
field["maximum_lag"],
field["maximum_lead"],
field["orig_maximum_endo_lag"],
field["orig_maximum_endo_lead"],
field["orig_maximum_exo_lag"],
field["orig_maximum_exo_lead"],
field["orig_maximum_exo_det_lag"],
field["orig_maximum_exo_det_lead"],
max(field["orig_maximum_lag"],
field["orig_maximum_lag_with_diffs_expanded"]),
field["orig_maximum_lead"]
)
function initialize_parameter!(params, field, symboltable)
s = symboltable[field["name"]]
k = s.orderintype
......@@ -135,10 +255,24 @@ function initialize_parameter!(params, field, symboltable)
end
function native_statement(field)
eval
# println("NATIVE: $field")
end
function initval(field)
function histval!(context, field, symboltable)
m = context.models[1]
histval = zeros(m.orig_maximum_lag, m.endogenous_nbr)
for v in field["vals"]
k = symboltable[v["name"]].orderintype
l = m.orig_maximum_lag - v["lag"]
histval[l, k] = Meta.parse(v["value"])
end
@show histval
context.work.histval = histval
@show context.work.histval
end
function initval!(field)
end
function shocks!(Sigma, field, symboltable)
......@@ -158,7 +292,8 @@ end
function set_stderr!(Sigma, stderr, symboltable)
for s in stderr
k = symboltable[s["name"]].orderintype
Sigma[k, k] = eval(Meta.parse(s["stderr"]))
x = eval(Meta.parse(s["stderr"]))
Sigma[k, k] = x*x
end
end
......@@ -207,23 +342,37 @@ function stoch_simul!(context, field)
model = context.models[1]
options = context.options
results = context.results.model_results[1]
work = context.work
options["stoch_simul"] = Dict()
copy!(options["stoch_simul"], field["options"])
compute_stoch_simul!(context)
compute_variance!(context)
x = results.linearrationalexpectations.g1
vx = view(x, :, 1:size(x, 2) - 1)
display_stoch_simul(vx', "Coefficients of approximate solution function", context)
if (periods = get(options["stoch_simul"], "periods", 0)) > 0
simulresults = Matrix{Float64}(undef, periods + 1, model.endogenous_nbr)
y0 = results.endogenous_steady_state
histval = work.histval
if size(histval, 1) == 0
y0 = results.endogenous_steady_state
else
y0 = view(histval, 1, :)
end
C = cholesky(model.Sigma_e)
x = randn(periods, model.exogenous_nbr)*C.U
x = vcat(zeros(1, model.exogenous_nbr),
randn(periods, model.exogenous_nbr)*C.U)
c = results.endogenous_steady_state
A = zeros(model.endogenous_nbr, model.endogenous_nbr)
B = zeros(model.endogenous_nbr, model.exogenous_nbr)
make_A_B!(A, B, model, results)
simul_first_order!(simulresults, y0, x, c, A, B, periods)
@show simulresults
if work.model_has_trend
simulresults .+= transpose(results.endogenous_trend[:,1]) .+ collect(0:size(simulresults, 1)-1) * transpose(results.endogenous_trend[:, 2])
end
first_period = get(options["stoch_simul"], "first_periods", 1)
endogenous_names = [Symbol(n) for n in get_endogenous_longname(context.symboltable)]
data = TimeDataFrame(simulresults, Periods.Undated, first_period, endogenous_names)
push!(results.simulations, Simulation("", "stoch_simul", options["stoch_simul"], data))
end
end
......@@ -319,3 +468,160 @@ function get_jacobian_at_steadystate!(work::Work, steadystate, exogenous, m::Mod
period)
end
function compute_variance!(context)
m = context.models[1]
results = context.results.model_results[1]
Σy = results.endogenous_variance
A = results.linearrationalexpectations.gs1
B1 = zeros(m.n_states, m.exogenous_nbr)
g1_1 = results.linearrationalexpectations.g1_1
g1_2 = results.linearrationalexpectations.g1_2
vr1 = view(g1_2, m.i_bkwrd_b, :)
B1 .= vr1
ws = LyapdWs(m.n_states)
Σ = zeros(m.n_states, m.n_states)
tmp = zeros(m.n_states, m.exogenous_nbr)
mul!(tmp, B1, m.Sigma_e)
B = zeros(m.n_states, m.n_states)
mul!(B, tmp, B1')
extended_lyapd!(Σ, A, B, ws)
vΣy = view(Σy, m.i_bkwrd_b, m.i_bkwrd_b)
vΣy .= Σ
n_non_states = m.endogenous_nbr - m.n_states
B2 = zeros(n_non_states, m.exogenous_nbr)
vr2 = view(results.linearrationalexpectations.g1_2, m.i_non_states, :)
B2 .= vr2
Σ = zeros(m.endogenous_nbr, n_non_states)
A2 = zeros(n_non_states, m.n_states)
vr3 = view(results.linearrationalexpectations.g1_1, m.i_non_states, :)
A2 .= vr3
tmp1 = zeros(m.endogenous_nbr, m.n_states)
mul!(tmp1, g1_1, vΣy)
mul!(Σ, tmp1, A2')
tmp2 = zeros(m.endogenous_nbr, m.exogenous_nbr)
mul!(tmp2, g1_2, m.Sigma_e)
mul!(Σ, tmp2, B2', 1.0, 1.0)
display(Σ)
vΣy = view(Σy, :, m.i_non_states)
vΣy .= Σ
vΣy = view(Σy, m.i_non_states, m.i_bkwrd_b)
= view(Σ, m.i_bkwrd_b, :)
vΣy .= '
end
function calib_smoother!(context, field, varobs, varobs_ids)
model = context.models[1]
options = context.options
results = context.results.model_results[1]
options["calib_smoother"] = Dict()
copy!(options["calib_smoother"], field["options"])
file = get(options["calib_smoother"], "datafile", "")
if (file = get(options["calib_smoother"], "datafile", "")) != ""
df = DataFrame(CSV.File(file))
if uppercase(names(df)[1]) in ["COLUMN1", "DATE", "DATES",
"PERIOD", "PERIODS"]
frequency = identify_period_frequency(uppercase(df[1, 1]))
else
frequency = Periods.Undated
firstperiod = 1
end
timedataframe = TimeDataFrame(df, frequency, firstperiod)
end
ny = length(varobs)
nobs = size(df, 1) - 1
start = 1
last = nobs
Y = Matrix{Float64}(undef, ny, nobs)
for (i, v) in enumerate(varobs)
Y[i, :] .= df[start + 1:last + 1, v] .- results.endogenous_steady_state[varobs_ids[i]]
end
statevar_ids = model.i_bkwrd_b
kalman_statevar_ids = collect(1:model.endogenous_nbr)
ns = length(kalman_statevar_ids)
np = model.exogenous_nbr
kws = KalmanSmootherWs{Float64, Int64}(ny, ns, model.exogenous_nbr, nobs)
c = zeros(ny)
k1 = findall(in(varobs_ids), kalman_statevar_ids)
k2 = findall(in(statevar_ids), kalman_statevar_ids)
Z = zeros(ny, ns)
for i in varobs_ids
Z[i, varobs_ids[1]] = 1.0
end
H = zeros(ny, ny)
d = zeros(ns)
T = zeros(ns, ns)
vg1 = view(context.results.model_results[1].linearrationalexpectations.g1_1, kalman_statevar_ids, :)
T[:, k2] .= vg1
R = zeros(ns, np)
vg2 = view(context.results.model_results[1].linearrationalexpectations.g1_2, kalman_statevar_ids, :)
R .= vg2
Q = model.Sigma_e
a0 = zeros(ns)
alphah = zeros(ns, nobs)
epsilonh = zeros(ny, nobs)
etah = zeros(np, nobs)
P = zeros(ns, ns)
vv = view(context.results.model_results[1].endogenous_variance, kalman_statevar_ids, kalman_statevar_ids)
P .= vv
Valpha = zeros(ns, ns, nobs)
Vepsilon = zeros(ny, ny, nobs)
Veta = zeros(np, np, nobs)
presample = 0
data_pattern = Vector{Vector{Integer}}(undef, 0)
for i = 1:nobs
push!(data_pattern, findall(Y[:, i] .!= NaN))
end
kalman_smoother!(Y, c, Z, H, d, T, R, Q, a0, P, alphah, epsilonh, etah,
Valpha, Vepsilon, Veta, start, last, presample,
kws, data_pattern)
alphah .+= results.endogenous_steady_state
results.smoother["alphah"] = alphah
Plots.plot(alphah, layout = 7)
end
function identify_period_frequency(period)::Periods.Frequency
if !isuppercase(period)
period = uppercase(period)
end
if 'Y' in period
frequency = Periods.Year
elseif 'A' in period
frequency = Periods.Year
elseif 'S' in period
frequency = Period.Semester
elseif 'H' in period
frequency = Periods.Semester
elseif 'Q' in period
frequency = Periods.Quarter
elseif 'M' in period
frequency = Periods.Month
elseif 'W' in period
frequency = Periods.Week
elseif 'D' in period
frequency = Periods.Day
elseif (cdash = count("-", period)) == 2
frequency = Periods.Day
elseif cdash == 1
frequency = Periods.Month
else
throw(ErrorException)
end
end
function deterministic_trends!(context, field)
model = context.models[1]
work = context.work
trend = context.results.model_results[1].endogenous_trend
symboltable = context.symboltable
work.model_has_trend = true
fill!(trend, 0.0)
for (key, value) in field["trends"]
trend[symboltable[key].orderintype,2] = work.params[symboltable[value].orderintype]
end