Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

nonlinear conditional forecasting #61

Merged
merged 11 commits into from
Dec 22, 2023
Merged

nonlinear conditional forecasting #61

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
36b7a36
nonlinear cond fcst; pruned sol no longer updating
thorek1 Dec 21, 2023
02d9bbc
update docstring
thorek1 Dec 21, 2023
0adb20e
fix cond fcst test
thorek1 Dec 21, 2023
7d17d6e
fix JET
thorek1 Dec 21, 2023
2369908
fix get_irf the 1000th time
thorek1 Dec 21, 2023
06b7d26
fix levels of cond fcst
thorek1 Dec 21, 2023
344a7cb
bundle get_relevant_steady_states
thorek1 Dec 21, 2023
6c059f5
fix plot_cond_fcst
thorek1 Dec 21, 2023
e281680
fix get_irf due to types
thorek1 Dec 22, 2023
93acdde
fix plotting
thorek1 Dec 22, 2023
513f110
tidy up
thorek1 Dec 22, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions docs/src/unfinished_docs/todo.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
- [ ] add technical details about SS solver, obc solver, and other algorithms
- [ ] fix translate dynare mod file from file written using write to dynare file (see test models)
- [ ] rm obc vars from get_SS
- [ ] write tests/docs for nonlinear obc and forecasting
- [ ] fix SS solver (failed for backus in guide)
- [ ] functions to reverse state_update (input: previous shock and current state, output previous state), find shocks corresponding to bringing one state to the next
- [ ] cover nested case: min(50,a+b+max(c,10))
Expand Down
128 changes: 100 additions & 28 deletions src/MacroModelling.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -303,6 +303,50 @@ function obc_objective_optim_fun(X::Vector{S}, grad::Vector{S}) where S
sum(abs2, X)
end



function match_conditions(res::Vector{S}, X::Vector{S}, jac::Matrix{S}, p) where S
Conditions, State_update, Shocks, Cond_var_idx, Free_shock_idx, State, Pruning, 𝒷 = p

if length(jac) > 0
jac .= 𝒜.jacobian(𝒷(), xx -> begin
Shocks[Free_shock_idx] .= xx

new_State = State_update(State, convert(typeof(xx), Shocks))

cond_vars = Pruning ? sum(new_State) : new_State

return abs2.(Conditions[Cond_var_idx] - cond_vars[Cond_var_idx])
end, X)[1]'
end

Shocks[Free_shock_idx] .= X

new_State = State_update(State, convert(typeof(X), Shocks))

cond_vars = Pruning ? sum(new_State) : new_State

res .= abs2.(Conditions[Cond_var_idx] - cond_vars[Cond_var_idx])
end



# function match_conditions(res::Vector{S}, X::Vector{S}, jac::Matrix{S}, p) where S
# conditions, state_update, shocks, cond_var_idx, free_shock_idx, state, 𝒷 = p

# if length(jac) > 0
# jac .= 𝒜.jacobian(𝒷(), xx -> begin
# shocks[free_shock_idx] .= xx
# return abs2.(conditions[cond_var_idx] - state_update(state, convert(typeof(xx), shocks))[cond_var_idx])
# end, X)[1]'
# end

# shocks[free_shock_idx] .= X

# res .= abs2.(conditions[cond_var_idx] - state_update(state, convert(typeof(X), shocks))[cond_var_idx])
# end


function set_up_obc_violation_function!(𝓂)
present_varss = collect(reduce(union,match_pattern.(get_symbols.(𝓂.dyn_equations),r"₍₀₎$")))

Expand Down Expand Up @@ -1088,6 +1132,30 @@ function parse_occasionally_binding_constraints(equations_block; max_obc_horizon
end



function get_relevant_steady_states(𝓂::ℳ, algorithm::Symbol)::Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}
full_NSSS = sort(union(𝓂.var,𝓂.aux,𝓂.exo_present))

full_NSSS[indexin(𝓂.aux,full_NSSS)] = map(x -> Symbol(replace(string(x), r"ᴸ⁽⁻?[⁰¹²³⁴⁵⁶⁷⁸⁹]+⁾" => "")), 𝓂.aux)

if any(x -> contains(string(x), "◖"), full_NSSS)
full_NSSS_decomposed = decompose_name.(full_NSSS)
full_NSSS = [length(a) > 1 ? string(a[1]) * "{" * join(a[2],"}{") * "}" * (a[end] isa Symbol ? string(a[end]) : "") : string(a[1]) for a in full_NSSS_decomposed]
end

relevant_SS = get_steady_state(𝓂, algorithm = algorithm, return_variables_only = true, derivatives = false)

reference_steady_state = [s ∈ 𝓂.exo_present ? 0 : relevant_SS(s) for s in full_NSSS]

relevant_NSSS = get_steady_state(𝓂, algorithm = :first_order, return_variables_only = true, derivatives = false)

NSSS = [s ∈ 𝓂.exo_present ? 0 : relevant_NSSS(s) for s in full_NSSS]

SSS_delta = NSSS - reference_steady_state

return reference_steady_state, NSSS, SSS_delta
end

# compatibility with SymPy
Max = max
Min = min
Expand Down Expand Up @@ -3264,10 +3332,7 @@ function solve!(𝓂::ℳ;
aug_state₁ = [pruned_states[1][𝓂.timings.past_not_future_and_mixed_idx]; 1; shock]
aug_state₂ = [pruned_states[2][𝓂.timings.past_not_future_and_mixed_idx]; 0; zero(shock)]

pruned_states[1] .= 𝐒₁ * aug_state₁
pruned_states[2] .= 𝐒₁ * aug_state₂ + 𝐒₂ * ℒ.kron(aug_state₁, aug_state₁) / 2

return pruned_states[1] + pruned_states[2] # strictly following Andreasen et al. (2018)
return [𝐒₁ * aug_state₁, 𝐒₁ * aug_state₂ + 𝐒₂ * ℒ.kron(aug_state₁, aug_state₁) / 2] # strictly following Andreasen et al. (2018)
end

if obc
Expand Down Expand Up @@ -3334,11 +3399,7 @@ function solve!(𝓂::ℳ;

kron_aug_state₁ = ℒ.kron(aug_state₁, aug_state₁)

pruned_states[1] .= 𝐒₁ * aug_state₁
pruned_states[2] .= 𝐒₁ * aug_state₂ + 𝐒₂ * kron_aug_state₁ / 2
pruned_states[3] .= 𝐒₁ * aug_state₃ + 𝐒₂ * ℒ.kron(aug_state₁̂, aug_state₂) + 𝐒₃ * ℒ.kron(kron_aug_state₁,aug_state₁) / 6

return pruned_states[1] + pruned_states[2] + pruned_states[3]
return [𝐒₁ * aug_state₁, 𝐒₁ * aug_state₂ + 𝐒₂ * kron_aug_state₁ / 2, 𝐒₁ * aug_state₃ + 𝐒₂ * ℒ.kron(aug_state₁̂, aug_state₂) + 𝐒₃ * ℒ.kron(kron_aug_state₁,aug_state₁) / 6]
end

if obc
Expand Down Expand Up @@ -4743,8 +4804,6 @@ function calculate_third_order_solution(∇₁::AbstractMatrix{<: Real}, #first
end




function irf(state_update::Function,
obc_state_update::Function,
initial_state::Union{Vector{Vector{Float64}},Vector{Float64}},
Expand Down Expand Up @@ -4939,22 +4998,30 @@ function irf(state_update::Function,

Y = zeros(T.nVars,periods,1)

Y[:,1,1] = state_update(initial_state,shock_history[:,1])
initial_state = state_update(initial_state,shock_history[:,1])

Y[:,1,1] = pruning ? sum(initial_state) : initial_state

for t in 1:periods-1
Y[:,t+1,1] = state_update(pruning ? initial_state : Y[:,t,1],shock_history[:,t+1])
initial_state = state_update(initial_state,shock_history[:,t+1])

Y[:,t+1,1] = pruning ? sum(initial_state) : initial_state
end

return KeyedArray(Y[var_idx,:,:] .+ level[var_idx]; Variables = axis1, Periods = 1:periods, Shocks = [:simulate])
elseif shocks == :none
Y = zeros(T.nVars,periods,1)

shck = T.nExo == 0 ? Vector{Float64}(undef, 0) : zeros(T.nExo)

Y[:,1,1] = state_update(initial_state,shck)

initial_state = state_update(initial_state, shck)

Y[:,1,1] = pruning ? sum(initial_state) : initial_state

for t in 1:periods-1
Y[:,t+1,1] = state_update(pruning ? initial_state : Y[:,t,1], shck)
initial_state = state_update(initial_state, shck)

Y[:,t+1,1] = pruning ? sum(initial_state) : initial_state
end

return KeyedArray(Y[var_idx,:,:] .+ level[var_idx]; Variables = axis1, Periods = 1:periods, Shocks = [:none])
Expand All @@ -4969,10 +5036,14 @@ function irf(state_update::Function,
shock_history[ii,1] = negative_shock ? -1 : 1
end

Y[:,1,i] = state_update(initial_state_copy, shock_history[:,1])
initial_state_copy = state_update(initial_state_copy, shock_history[:,1])

Y[:,1,i] = pruning ? sum(initial_state_copy) : initial_state_copy

for t in 1:periods-1
Y[:,t+1,i] = state_update(pruning ? initial_state_copy : Y[:,t,i], shock_history[:,t+1])
initial_state_copy = state_update(initial_state_copy, shock_history[:,t+1])

Y[:,t+1,i] = pruning ? sum(initial_state_copy) : initial_state_copy
end
end

Expand Down Expand Up @@ -5043,11 +5114,7 @@ function girf(state_update::Function,
initial_state_copy² = deepcopy(initial_state_copy)

for i in 1:warmup_periods
if pruning
state_update(initial_state_copy², randn(T.nExo))
else
initial_state_copy² = state_update(initial_state_copy², randn(T.nExo))
end
initial_state_copy² = state_update(initial_state_copy², randn(T.nExo))
end

Y₁ = zeros(T.nVars, periods + 1)
Expand All @@ -5061,11 +5128,13 @@ function girf(state_update::Function,
end

if pruning
Y₁[:,1] = state_update(initial_state_copy², baseline_noise)
Y₂[:,1] = state_update(initial_state_copy², baseline_noise)
initial_state_copy² = state_update(initial_state_copy², baseline_noise)

initial_state₁ = deepcopy(initial_state_copy²)
initial_state₂ = deepcopy(initial_state_copy²)

Y₁[:,1] = initial_state_copy² |> sum
Y₂[:,1] = initial_state_copy² |> sum
else
Y₁[:,1] = state_update(initial_state_copy², baseline_noise)
Y₂[:,1] = state_update(initial_state_copy², baseline_noise)
Expand All @@ -5075,8 +5144,11 @@ function girf(state_update::Function,
baseline_noise = randn(T.nExo)

if pruning
Y₁[:,t+1] = state_update(initial_state₁, baseline_noise)
Y₂[:,t+1] = state_update(initial_state₂, baseline_noise + shock_history[:,t])
initial_state₁ = state_update(initial_state₁, baseline_noise)
initial_state₂ = state_update(initial_state₂, baseline_noise + shock_history[:,t])

Y₁[:,t+1] = initial_state₁ |> sum
Y₂[:,t+1] = initial_state₂ |> sum
else
Y₁[:,t+1] = state_update(Y₁[:,t],baseline_noise)
Y₂[:,t+1] = state_update(Y₂[:,t],baseline_noise + shock_history[:,t])
Expand Down Expand Up @@ -6159,7 +6231,7 @@ function filter_and_smooth(𝓂::ℳ, data_in_deviations::AbstractArray{Float64}
PCiF = P[:, :, t] * C' * iF[:, :, t]
L[:, :, t] .= A - A * PCiF * C
P[:, :, t+1].= A * P[:, :, t] * L[:, :, t]' + 𝐁
σ[:, t] .= sqrt.(abs.(ℒ.diag(P[:, :, t+1]))) # small numerica errors in this computation
σ[:, t] .= sqrt.(abs.(ℒ.diag(P[:, :, t+1]))) # small numerical errors in this computation
μ[:, t+1] .= A * (μ[:, t] + PCiF * v[:, t])
ϵ[:, t] .= B' * C' * iF[:, :, t] * v[:, t]
end
Expand Down
2 changes: 1 addition & 1 deletion src/common_docstrings.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ const NEGATIVE_SHOCK = "`negative_shock` [Default: `false`, Type: `Bool`]: calcu
const GENERALISED_IRF = "`generalised_irf` [Default: `false`, Type: `Bool`]: calculate generalised IRFs. Relevant for nonlinear solutions. Reference steady state for deviations is the stochastic steady state."
const INITIAL_STATE = "`initial_state` [Default: `[0.0]`, Type: `Vector{Float64}`]: provide state (in levels, not deviations) from which to start IRFs. Relevant for normal IRFs. The state includes all variables as well as exogenous variables in leads or lags if present."
const ALGORITHM = "`algorithm` [Default: `:first_order`, Type: `Symbol`]: algorithm to solve for the dynamics of the model."
const LEVELS = "`levels` [Default: `false`, Type: `Bool`]: return levels or absolute deviations from steady state."
const LEVELS = "`levels` [Default: `false`, Type: `Bool`]: return levels or absolute deviations from steady state corresponding to the solution algorithm (e.g. stochastic steady state for higher order solution algorithms)."
const CONDITIONS = "`conditions` [Type: `Union{Matrix{Union{Nothing,Float64}}, SparseMatrixCSC{Float64}, KeyedArray{Union{Nothing,Float64}}, KeyedArray{Float64}}`]: conditions for which to find the corresponding shocks. The input can have multiple formats, but for all types of entries the first dimension corresponds to the number of variables and the second dimension to the number of periods. The conditions can be specified using a matrix of type `Matrix{Union{Nothing,Float64}}`. In this case the conditions are matrix elements of type `Float64` and all remaining (free) entries are `nothing`. You can also use a `SparseMatrixCSC{Float64}` as input. In this case only non-zero elements are taken as conditions. Note that you cannot condition variables to be zero using a `SparseMatrixCSC{Float64}` as input (use other input formats to do so). Another possibility to input conditions is by using a `KeyedArray`. You can use a `KeyedArray{Union{Nothing,Float64}}` where, similar to `Matrix{Union{Nothing,Float64}}`, all entries of type `Float64` are recognised as conditions and all other entries have to be `nothing`. Furthermore, you can specify in the primary axis a subset of variables (of type `Symbol` or `String`) for which you specify conditions and all other variables are considered free. The same goes for the case when you use `KeyedArray{Float64}}` as input, whereas in this case the conditions for the specified variables bind for all periods specified in the `KeyedArray`, because there are no `nothing` entries permitted with this type."
const SHOCK_CONDITIONS = "`shocks` [Default: `nothing`, Type: `Union{Matrix{Union{Nothing,Float64}}, SparseMatrixCSC{Float64}, KeyedArray{Union{Nothing,Float64}}, KeyedArray{Float64}, Nothing} = nothing`]: known values of shocks. This entry allows the user to include certain shock values. By entering restrictions on the shock sin this way the problem to match the conditions on endogenous variables is restricted to the remaining free shocks in the repective period. The input can have multiple formats, but for all types of entries the first dimension corresponds to the number of shocks and the second dimension to the number of periods. The shocks can be specified using a matrix of type `Matrix{Union{Nothing,Float64}}`. In this case the shocks are matrix elements of type `Float64` and all remaining (free) entries are `nothing`. You can also use a `SparseMatrixCSC{Float64}` as input. In this case only non-zero elements are taken as certain shock values. Note that you cannot condition shocks to be zero using a `SparseMatrixCSC{Float64}` as input (use other input formats to do so). Another possibility to input known shocks is by using a `KeyedArray`. You can use a `KeyedArray{Union{Nothing,Float64}}` where, similar to `Matrix{Union{Nothing,Float64}}`, all entries of type `Float64` are recognised as known shocks and all other entries have to be `nothing`. Furthermore, you can specify in the primary axis a subset of shocks (of type `Symbol` or `String`) for which you specify values and all other shocks are considered free. The same goes for the case when you use `KeyedArray{Float64}}` as input, whereas in this case the values for the specified shocks bind for all periods specified in the `KeyedArray`, because there are no `nothing` entries permitted with this type."
const PARAMETER_DERIVATIVES = "`parameter_derivatives` [Default: :all]: parameters for which to calculate partial derivatives. Inputs can be a parameter name passed on as either a `Symbol` or `String` (e.g. `:alpha`, or \"alpha\"), or `Tuple`, `Matrix` or `Vector` of `String` or `Symbol`. `:all` will include all parameters."
Expand Down
Loading
Loading