Add balanced growth path support for first-order models

[ylevch]

Hi Thore,

I wanted to be able to work with models with levels and BGP, so I took a stab at it.
This PR adds first-order support for such models with auto-detection in the equation parser.

I am a relative novice to open source development/contribution process and learn many things as I go (with Claude's help). Feedback would be appreciated. As a reference, I used IRIS-Toolbox implementation, so the main idea is that every variable gets a level and an additive growth unknown x^G. Since we get 2N unknowns, each steady-state equation is evaluated at two time origins (t, t+1) to double the equation count.

This PR includes:

• Auto-detection of non-stationary levels in the parser; build the 2N augmented (level + growth) steady-state system; anchor the free level of each independent cointegrated trend to 0 via the steady-state Jacobian's rank.
• ∂SS/∂param inverts a singular augmented Jacobian (free trend levels give zero columns; the two-time-point construction gives duplicate rows), fixed with the minimum-norm least-squares solution (pinv).
• get_solution/get_irf work end-to-end for levels models.
• Deterministic BGP drift x_t = anchor + xᴳ·t is added to trending variables in get_irf/simulate when levels = true. Unconditional variances of trending levels are NaN.
• x[ss] = expr lets one to pin a trend's level explicitly (otherwise automatic anchoring is to 0).
• I also included toy tests and compared steady state outputs of a simple QPM using this implementation vs IRIS-Toolbox: levels, growth rates and cointegration relationships matched.

Best,
Lazar

[thorek1]

Hi Lazar,

This is really cool!
Thank you for spending time on this. Implementing BGP handling has been on my todo list since a long time and I was thinking of going down this route.

I had a brief look at your PR and my only comment at this stage would be to include a more complex model, ideally a benchmark test case for this kind of problem.

I will go through the codes in more detail in the coming weeks.

Best,

Thore

[codecov-commenter]

Codecov Report

❌ Patch coverage is 92.30769% with 13 lines in your changes missing coverage. Please review.
✅ Project coverage is 32.18%. Comparing base (f2274a0) to head (70bd494).

Files with missing lines	Patch %	Lines
ext/ForwardDiffExt.jl	0.00%	8 Missing ⚠️
src/rrules.jl	50.00%	3 Missing ⚠️
src/parser/equation_processing.jl	96.42%	1 Missing ⚠️
src/parser/model_setup.jl	97.61%	1 Missing ⚠️

❗ There is a different number of reports uploaded between BASE (f2274a0) and HEAD (70bd494). Click for more details.

HEAD has 31 uploads less than BASE

Flag	BASE (f2274a0)	HEAD (70bd494)
	32	1

Additional details and impacted files

@@             Coverage Diff             @@
##             main     #304       +/-   ##
===========================================
- Coverage   88.05%   32.18%   -55.87%     
===========================================
  Files          35       35               
  Lines       30243    29372      -871     
===========================================
- Hits        26629     9453    -17176     
- Misses       3614    19919    +16305     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[thorek1]

I will collect some more thoughts here as I go along looking at the PR.

Suggestions:

• covariance of BGP models refers to the covariance along the path (finite values instead of the currently infinite ones)
• Show the growth rates in a separate column after the steady state values in SS calls
• make it work with models with expectations. e.g.

using MacroModelling

@model RBC_baseline begin
	c[0] ^ (-σ) = β * c[1] ^ (-σ) * (α * z[1] * (k[0] / l[1]) ^ (α - 1) + 1 - δ)

	ψ * c[0] ^ σ / (1 - l[0]) = w[0]

	k[0] = (1 - δ) * k[-1] + i[0]

	y[0] = c[0] + i[0] + g[0]

	y[0] = z[0] * k[-1] ^ α * l[0] ^ (1 - α)

	w[0] = y[0] * (1 - α) / l[0]

	r[0] = y[0] * α * 4 / k[-1]

	z[0] = (1 - ρᶻ) + ρᶻ * z[-1] + σᶻ * ϵᶻ[x]

	g[0] = (1 - ρᵍ) * ḡ + ρᵍ * g[-1] + σᵍ * ϵᵍ[x]

end


@parameters RBC_baseline begin
	σᶻ = 0.066

	σᵍ = .104

	σ = 1

	α = 1/3

	i_y = 0.25

	k_y = 10.4

	ρᶻ = 0.97

	ρᵍ = 1.01 # unit root

	g_y = 0.2038

	ḡ | ḡ = g_y * y[ss]

    δ = i_y / k_y

    β = 1 / (α / k_y + (1 - δ))

	ψ | l[ss] = 1/3
end

get_solution(RBC_baseline)