diff --git a/README.md b/README.md
index e94b4cc2d..8974a9264 100644
--- a/README.md
+++ b/README.md
@@ -141,7 +141,7 @@ The package contains the following models in the `models` folder:
 **Occasionally binding constraints**|yes|yes|yes|yes|yes||||yes|||yes||
 **Global solution**||||yes|yes|||||||yes||
 **Estimation**|yes|yes|yes|||yes||yes|yes|yes|yes|||
-**Balanced growth path**||yes|yes||||yes|yes|yes|yes|||||
+**Balanced growth path**|yes|yes|yes||||yes|yes|yes|yes|||||
 **Model input**|macro (julia)|text file|text file|text file|text file|macro (julia)|module (julia)|text file|text file|text file|text file|text file|text file|
 **Timing convention**|end-of-period|end-of-period||end-of-period|start-of-period|start-of-period|end-of-period|end-of-period|end-of-period|end-of-period|end-of-period|start-of-period|start-of-period|
 
diff --git a/docs/src/how-to/balanced_growth.md b/docs/src/how-to/balanced_growth.md
new file mode 100644
index 000000000..fb2f19bb7
--- /dev/null
+++ b/docs/src/how-to/balanced_growth.md
@@ -0,0 +1,184 @@
+# Balanced Growth Path Handling
+
+Many DSGE models feature non-stationary variables that grow along a balanced growth path (BGP). Common examples include output, consumption, capital, and wages in models with technological progress. `MacroModelling.jl` provides automatic handling of such models through the `auto_detrend`, `deflator`, and `trend_var` options.
+
+## Overview
+
+In a model with balanced growth, certain variables grow at constant rates over time. For example, with labor-augmenting technological progress at rate `γ`, output, consumption, and capital all grow at rate `γ` in the long run.
+
+To solve such models using perturbation methods, we need to transform the non-stationary variables into stationary ones by "detrending" - dividing by the level of technology (or another appropriate trend variable).
+
+`MacroModelling.jl` automates this process with two approaches:
+
+### Automatic Detection (Recommended)
+
+Use `auto_detrend = true` to let the package automatically:
+1. Detect trend variables from equation patterns (e.g., `A[0] = γ * A[-1]`)
+2. Identify which variables should be deflated by each trend
+3. Apply the detrending transformation
+4. Solve for the detrended steady state
+
+### Manual Specification
+
+Alternatively, you can manually specify:
+1. Which variables should be detrended using `deflator` in `@model`
+2. Trend variable growth factors using `trend_var` in `@parameters`
+
+## Automatic Detrending
+
+The simplest way to handle balanced growth is to use automatic detection:
+
+```julia
+@model RBC_growth auto_detrend = true begin
+    # Technology grows at rate γ - automatically detected as trend variable
+    A[0] = γ * A[-1]
+    
+    # Variables multiplied by A will be auto-detected as trending
+    y[0] = A[0] * k[-1]^α
+    c[0] + k[0] = y[0] + (1-δ)*k[-1]
+    1/c[0] = β * (1/c[1]) * (α * y[1]/k[0] + 1-δ)
+end
+
+@parameters RBC_growth begin
+    γ = 1.02    # 2% growth rate
+    α = 0.33
+    β = 0.99
+    δ = 0.025
+end
+```
+
+The package will:
+1. Detect `A` as a trend variable (from `A[0] = γ * A[-1]`)
+2. Identify `y`, `c`, `k` as variables that should be divided by `A`
+3. Transform the equations to work with detrended (stationary) variables
+4. Solve for the balanced growth path steady state
+
+## Manual Specification
+
+### Specifying Deflators
+
+Use the `deflator` option in `@model` to specify which variables are non-stationary and what their deflator (trend variable) is:
+
+```julia
+@model RBC_growth deflator = Dict(:y => :A, :c => :A, :k => :A) begin
+    # A is the trend variable (e.g., technology level)
+    A[0] = γ * A[-1]
+    
+    # Write equations in terms of detrended variables (y/A, c/A, k/A)
+    # The package transforms y[0] → y[0] * A[0] internally
+    y[0] = k[-1]^α
+    c[0] + k[0] = y[0] + (1-δ)*k[-1]
+    1/c[0] = β * (1/c[1]) * (α * y[1]/k[0] + 1-δ)
+end
+```
+
+The `deflator` option takes a `Dict{Symbol, Symbol}` where:
+- Keys are the non-stationary variables to be detrended
+- Values are the trend variables to use as deflators
+
+### Specifying Trend Variable Growth Factors
+
+Use the `trend_var` option in `@parameters` to document the growth factors of trend variables:
+
+```julia
+@parameters RBC_growth trend_var = Dict(:A => :γ) begin
+    γ = 1.02    # 2% growth rate
+    α = 0.33
+    β = 0.99
+    δ = 0.025
+end
+```
+
+## How It Works
+
+### Automatic Detection
+
+When `auto_detrend = true`, the package:
+
+1. **Identifies trend variables**: Scans equations for patterns like `X[0] = g * X[-1]` or `X[0] = X[-1] * g` where `g` is a parameter or expression. These indicate unit root / deterministic trend processes.
+
+2. **Identifies trending variables**: Finds variables that appear in the same equations as trend variables (heuristic approach).
+
+3. **Applies detrending**: For each trending variable `v` with trend `T`, transforms `v[t] → v[t] / T[t]`.
+
+### Manual Specification
+
+When you specify a deflator for a variable, `MacroModelling.jl` automatically transforms the equations. For each variable `v` with deflator `d`:
+
+- `v[0]` is transformed to `(v[0] * d[0])`
+- `v[-1]` is transformed to `(v[-1] * d[-1])`
+- `v[1]` is transformed to `(v[1] * d[1])`
+- `v[ss]` is transformed to `(v[ss] * d[ss])`
+
+This means you write your model in terms of detrended variables, and the package automatically "re-trends" them to recover the original (level) equations.
+
+## Checking Balanced Growth Configuration
+
+You can inspect the balanced growth path configuration using these functions:
+
+```julia
+# Check if a model has balanced growth handling enabled
+has_balanced_growth(model)
+
+# Get detailed information about the balanced growth configuration
+info = get_balanced_growth_path_info(model)
+# Returns: (has_balanced_growth, trend_vars, deflators, detrended_vars)
+```
+
+## Example: RBC Model with Technological Progress
+
+Here's a complete example of an RBC model with exogenous labor-augmenting technological progress using automatic detection:
+
+```julia
+using MacroModelling
+
+# Define the model with automatic detrending
+@model RBC_BGP auto_detrend = true begin
+    # Technology grows at rate γ
+    A[0] = γ * A[-1]
+    
+    # Production function
+    y[0] = A[0] * k[-1]^α
+    
+    # Resource constraint
+    c[0] + k[0] = y[0] + (1-δ)*k[-1]
+    
+    # Euler equation
+    1/c[0] = β * (1/c[1]) * (α * y[1]/k[0] + 1-δ)
+    
+    # Productivity shock (stationary)
+    z[0] = ρ * z[-1] + σ * eps[x]
+end
+
+# Define parameters
+@parameters RBC_BGP begin
+    γ = 1.005   # 0.5% quarterly growth
+    α = 0.33
+    β = 0.99
+    δ = 0.025
+    ρ = 0.9
+    σ = 0.01
+end
+
+# The model can now be solved and analyzed
+get_steady_state(RBC_BGP)
+```
+
+## Notes and Best Practices
+
+1. **Trend Variable Equation**: Always include an equation defining the evolution of the trend variable (e.g., `A[0] = γ * A[-1]`).
+
+2. **Consistent Detrending**: Make sure all variables that grow at the same rate use the same deflator.
+
+3. **Steady State**: The detrended variables should have well-defined steady states. If steady state solving fails, consider providing initial guesses.
+
+4. **Automatic vs Manual**: Use `auto_detrend = true` for simple cases. For complex models with multiple trends or specific detrending requirements, manual specification with `deflator` provides more control.
+
+5. **Comparison with Other Packages**: This feature is similar to Dynare's `trend_var` and `deflator` options, though the syntax differs slightly.
+
+## See Also
+
+- [`@model`](@ref) - Main model definition macro
+- [`@parameters`](@ref) - Parameter definition macro
+- [`get_balanced_growth_path_info`](@ref) - Get balanced growth configuration
+- [`has_balanced_growth`](@ref) - Check if balanced growth is enabled
diff --git a/docs/src/index.md b/docs/src/index.md
index 66f230727..d8c1ce484 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -63,7 +63,7 @@ The package contains the following models in the `models` folder:
 |**Occasionally binding constraints**|yes|yes|yes|yes|yes||||yes|||yes||
 |**Global solution**||||yes|yes|||||||yes||
 |**Estimation**|yes|yes|yes|||yes||yes|yes|yes|yes|||
-|**Balanced growth path**||yes|yes||||yes|yes|yes|yes||||
+|**Balanced growth path**|yes|yes|yes||||yes|yes|yes|yes||||
 |**Model input**|macro (julia)|text file|text file|text file|text file|macro (julia)|module (julia)|text file|text file|text file|text file|text file|text file|
 |**Timing convention**|end-of-period|end-of-period||end-of-period|start-of-period|start-of-period|end-of-period|end-of-period|end-of-period|end-of-period|end-of-period|start-of-period|start-of-period|
 
diff --git a/docs/src/unfinished_docs/todo.md b/docs/src/unfinished_docs/todo.md
index ceecd7285..d104006db 100644
--- a/docs/src/unfinished_docs/todo.md
+++ b/docs/src/unfinished_docs/todo.md
@@ -109,7 +109,7 @@
 - [ ] check warnings, errors throughout. check suppress not interfering with pigeons
 - [ ] 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))
-- [ ] add balanced growth path handling
+- [x] add balanced growth path handling
 - [ ] autocorr and covariance with derivatives. return 3d array
 - [ ] add pydsge and econpizza to overview
 - [ ] add for loop parser in @parameters
diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl
index b47e87f69..365794812 100644
--- a/src/MacroModelling.jl
+++ b/src/MacroModelling.jl
@@ -207,6 +207,10 @@ export translate_mod_file, translate_dynare_file, import_model, import_dynare
 export write_mod_file, write_dynare_file, write_to_dynare_file, write_to_dynare, export_dynare, export_to_dynare, export_mod_file, export_model
 
 export get_equations, get_steady_state_equations, get_dynamic_equations, get_calibration_equations, get_parameters, get_calibrated_parameters, get_parameters_in_equations, get_parameters_defined_by_parameters, get_parameters_defining_parameters, get_calibration_equation_parameters, get_variables, get_nonnegativity_auxiliary_variables, get_dynamic_auxiliary_variables, get_shocks, get_state_variables, get_jump_variables, get_missing_parameters, has_missing_parameters
+
+# Balanced growth path
+export has_balanced_growth, get_balanced_growth_info, get_trend_variables, get_variable_degree
+
 # Internal
 export irf, girf
 
diff --git a/src/inspect.jl b/src/inspect.jl
index 1855ba2c6..6c2825190 100644
--- a/src/inspect.jl
+++ b/src/inspect.jl
@@ -943,4 +943,113 @@ function get_jump_variables(𝓂::ℳ)::Vector{String}
     𝓂.timings.future_not_past_and_mixed |> collect |> sort .|> x -> replace.(string.(x), "◖" => "{", "◗" => "}")
 end
 
+
+"""
+$(SIGNATURES)
+Check if a model has balanced growth path structure (trend variables detected).
+
+# Arguments
+- `𝓂`: A model object
+
+# Returns
+- `Bool`: `true` if trend variables were detected, `false` otherwise.
+
+# Examples
+```julia
+using MacroModelling
+
+@model RBC_growth begin
+    A[0] = γ * A[-1]  # Trend variable
+    y[0] = A[0] * k[-1]^α
+    # ...
+end
+
+@parameters RBC_growth begin
+    γ = 1.02
+    α = 0.33
+end
+
+has_balanced_growth(RBC_growth)  # returns true
+```
+"""
+function has_balanced_growth(𝓂::ℳ)::Bool
+    !isempty(𝓂.balanced_growth.trend_variables)
+end
+
+
+"""
+$(SIGNATURES)
+Get information about the balanced growth path structure of the model.
+
+# Arguments
+- `𝓂`: A model object
+
+# Returns
+- `NamedTuple` with fields:
+  - `trend_variables`: Dict mapping trend variable names to their growth rate parameters
+  - `variable_degrees`: Dict mapping variable names to their homogeneity degrees
+  - `growth_parameters`: Set of parameters representing growth rates
+
+# Examples
+```julia
+using MacroModelling
+
+@model RBC_growth begin
+    A[0] = γ * A[-1]
+    y[0] = A[0] * k[-1]^α
+    c[0] + k[0] = y[0] + (1-δ)*k[-1]
+    1/c[0] = β * (1/c[1]) * (α * y[1]/k[0] + 1-δ)
+end
+
+@parameters RBC_growth begin
+    γ = 1.02
+    α = 0.33
+    β = 0.99
+    δ = 0.025
+end
+
+info = get_balanced_growth_info(RBC_growth)
+# info.trend_variables  # Dict(:A => :γ)
+# info.variable_degrees # Dict(:A => 1.0, :y => 1.0, :c => 1.0, :k => 1.0)
+```
+"""
+function get_balanced_growth_info(𝓂::ℳ)
+    (
+        trend_variables = 𝓂.balanced_growth.trend_variables,
+        variable_degrees = 𝓂.balanced_growth.variable_degrees,
+        growth_parameters = 𝓂.balanced_growth.growth_rate_parameters
+    )
+end
+
+
+"""
+$(SIGNATURES)
+Get the names of trend variables in the model.
+
+# Arguments
+- `𝓂`: A model object
+
+# Returns
+- `Vector{Symbol}`: Names of variables that have been identified as trends (unit roots).
+"""
+function get_trend_variables(𝓂::ℳ)::Vector{Symbol}
+    collect(keys(𝓂.balanced_growth.trend_variables)) |> sort
+end
+
+
+"""
+$(SIGNATURES)
+Get the homogeneity degree of a specific variable.
+
+# Arguments
+- `𝓂`: A model object
+- `var::Symbol`: The variable name
+
+# Returns
+- `Float64`: The homogeneity degree (0.0 for stationary, 1.0 for growing at trend rate)
+"""
+function get_variable_degree(𝓂::ℳ, var::Symbol)::Float64
+    get(𝓂.balanced_growth.variable_degrees, var, 0.0)
+end
+
 end # dispatch_doctor
\ No newline at end of file
diff --git a/src/macros.jl b/src/macros.jl
index c7c0f8ed1..83d5eccf1 100644
--- a/src/macros.jl
+++ b/src/macros.jl
@@ -1,6 +1,678 @@
 const all_available_algorithms = [:first_order, :second_order, :pruned_second_order, :third_order, :pruned_third_order]
 
 
+# ============================================================================
+# Balanced Growth Path Functions
+# ============================================================================
+
+"""
+    extract_variables_from_expr(expr)
+
+Extract all variable symbols from an expression that appear in `ref` expressions 
+(e.g., `var[0]`, `var[-1]`). Returns a Set{Symbol} of variable names.
+"""
+function extract_variables_from_expr(expr)
+    vars = Set{Symbol}()
+    postwalk(x -> begin
+        if x isa Expr && x.head == :ref
+            var_name = x.args[1]
+            if var_name isa Symbol
+                push!(vars, var_name)
+            end
+        end
+        x
+    end, expr)
+    return vars
+end
+
+
+"""
+    detect_trend_variables(model_ex::Expr)
+
+Detect trend (unit root) variables from model equations.
+
+A variable X is identified as a trend variable if it appears in an equation of the form:
+- `X[0] = γ * X[-1]` (deterministic trend with growth rate γ)
+- `X[0] = X[-1] * γ` (same pattern, different order)
+
+This pattern indicates X has a unit root and grows at rate γ over time.
+
+Returns a Dict mapping trend variable symbols to their growth rate parameter symbols.
+"""
+function detect_trend_variables(model_ex::Expr)
+    trend_vars = Dict{Symbol, Symbol}()
+    
+    for arg in model_ex.args
+        if !isa(arg, Expr) || arg.head != :(=)
+            continue
+        end
+        
+        lhs, rhs = arg.args[1], arg.args[2]
+        
+        # Check if lhs is var[0]
+        if !(lhs isa Expr && lhs.head == :ref && lhs.args[2] == 0)
+            continue
+        end
+        
+        var_name = lhs.args[1]
+        if !(var_name isa Symbol)
+            continue
+        end
+        
+        # Check if rhs is growth_factor * var[-1] or var[-1] * growth_factor
+        if rhs isa Expr && rhs.head == :call && rhs.args[1] == :*
+            term1 = rhs.args[2]
+            term2 = length(rhs.args) >= 3 ? rhs.args[3] : nothing
+            
+            if isnothing(term2)
+                continue
+            end
+            
+            # Pattern: var[-1] * growth_factor
+            if term1 isa Expr && term1.head == :ref && 
+               term1.args[1] == var_name && term1.args[2] == -1 &&
+               term2 isa Symbol
+                trend_vars[var_name] = term2
+            # Pattern: growth_factor * var[-1]
+            elseif term2 isa Expr && term2.head == :ref && 
+                   term2.args[1] == var_name && term2.args[2] == -1 &&
+                   term1 isa Symbol
+                trend_vars[var_name] = term1
+            end
+        end
+    end
+    
+    return trend_vars
+end
+
+
+"""
+    compute_homogeneity_degree(expr, var_degrees::Dict{Symbol, Float64}, 
+                               params::Set{Symbol}, trend_vars::Set{Symbol})
+
+Recursively compute the homogeneity degree of an expression.
+
+Homogeneity rules:
+- Variables have their assigned degree from `var_degrees`
+- Trend variables have degree 1.0
+- Parameters and constants have degree 0.0
+- For `a * b`: degree(a) + degree(b)
+- For `a / b`: degree(a) - degree(b)  
+- For `a^n` where n is numeric: n * degree(a)
+- For `a^p` where p is a parameter: degree(a) (assumes parameter is O(1))
+- For `a + b` or `a - b`: both must have same degree (returns that degree)
+- For functions like `exp`, `log`: argument must have degree 0
+
+Returns (degree::Float64, valid::Bool) where valid indicates if homogeneity holds.
+"""
+function compute_homogeneity_degree(expr, var_degrees::Dict{Symbol, Float64}, 
+                                    params::Set{Symbol}, trend_vars::Set{Symbol})
+    if expr isa Number
+        return (0.0, true)
+    elseif expr isa Symbol
+        if expr ∈ trend_vars
+            return (1.0, true)
+        elseif haskey(var_degrees, expr)
+            return (var_degrees[expr], true)
+        else
+            return (0.0, true)  # Parameter or constant
+        end
+    elseif expr isa Expr
+        if expr.head == :ref
+            var_name = expr.args[1]
+            if var_name isa Symbol
+                if var_name ∈ trend_vars
+                    return (1.0, true)
+                elseif haskey(var_degrees, var_name)
+                    return (var_degrees[var_name], true)
+                end
+            end
+            return (0.0, true)
+        elseif expr.head == :call
+            op = expr.args[1]
+            
+            if op == :*
+                total = 0.0
+                for i in 2:length(expr.args)
+                    d, valid = compute_homogeneity_degree(expr.args[i], var_degrees, params, trend_vars)
+                    if !valid
+                        return (0.0, false)
+                    end
+                    total += d
+                end
+                return (total, true)
+                
+            elseif op == :/
+                if length(expr.args) >= 3
+                    d1, v1 = compute_homogeneity_degree(expr.args[2], var_degrees, params, trend_vars)
+                    d2, v2 = compute_homogeneity_degree(expr.args[3], var_degrees, params, trend_vars)
+                    if !v1 || !v2
+                        return (0.0, false)
+                    end
+                    return (d1 - d2, true)
+                end
+                
+            elseif op == :^
+                base_d, base_v = compute_homogeneity_degree(expr.args[2], var_degrees, params, trend_vars)
+                if !base_v
+                    return (0.0, false)
+                end
+                exp_val = expr.args[3]
+                if exp_val isa Number
+                    return (base_d * exp_val, true)
+                elseif exp_val isa Symbol && exp_val ∈ params
+                    # For x^α where α is a parameter, we assume the result grows at base rate
+                    return (base_d, true)
+                else
+                    # Complex exponent - assume degree 0 for safety
+                    return (0.0, base_d == 0.0)
+                end
+                
+            elseif op in (:+, :-)
+                if length(expr.args) < 2
+                    return (0.0, true)
+                end
+                d1, v1 = compute_homogeneity_degree(expr.args[2], var_degrees, params, trend_vars)
+                if !v1
+                    return (0.0, false)
+                end
+                for i in 3:length(expr.args)
+                    di, vi = compute_homogeneity_degree(expr.args[i], var_degrees, params, trend_vars)
+                    if !vi || abs(di - d1) > 1e-10
+                        return (0.0, false)  # Mismatched degrees
+                    end
+                end
+                return (d1, true)
+                
+            elseif op in (:exp, :log, :sin, :cos, :tan, :sqrt)
+                d, v = compute_homogeneity_degree(expr.args[2], var_degrees, params, trend_vars)
+                # These functions require degree-0 arguments
+                return (0.0, v && abs(d) < 1e-10)
+            end
+        end
+    end
+    return (0.0, true)
+end
+
+
+"""
+    solve_homogeneity_constraints(model_ex::Expr, trend_vars::Dict{Symbol, Symbol}, 
+                                  all_vars::Set{Symbol}, params::Set{Symbol})
+
+Solve for the homogeneity degrees of all variables using iterative constraint propagation.
+
+The algorithm:
+1. Assign degree 1 to trend variables
+2. For each equation with form `LHS = RHS`, extract degree constraints
+3. Propagate degrees: if v appears in an expression with known-degree terms, infer v's degree
+4. Iterate until convergence
+
+For balanced growth models, trending variables (degree 1) include:
+- Output, consumption, capital when they grow with productivity
+- Any variable that appears multiplied by the trend or other trending variables
+
+Returns Dict{Symbol, Float64} mapping each variable to its homogeneity degree.
+"""
+function solve_homogeneity_constraints(model_ex::Expr, trend_vars::Dict{Symbol, Symbol},
+                                       all_vars::Set{Symbol}, params::Set{Symbol})
+    # Initialize degrees: trend vars have degree 1, others unknown (-1)
+    degrees = Dict{Symbol, Float64}()
+    trend_var_set = Set(keys(trend_vars))
+    
+    for v in trend_var_set
+        degrees[v] = 1.0
+    end
+    
+    # Non-trend variables start unknown
+    non_trend_vars = setdiff(all_vars, trend_var_set)
+    for v in non_trend_vars
+        degrees[v] = -1.0  # Unknown
+    end
+    
+    # Skip the trend equation itself when inferring degrees
+    trend_equations = Set{Int}()
+    for (i, arg) in enumerate(model_ex.args)
+        if isa(arg, Expr) && arg.head == :(=)
+            lhs = arg.args[1]
+            if lhs isa Expr && lhs.head == :ref
+                var_name = lhs.args[1]
+                if var_name isa Symbol && var_name ∈ trend_var_set
+                    push!(trend_equations, i)
+                end
+            end
+        end
+    end
+    
+    # Iteratively propagate degrees
+    max_iterations = 10
+    for _ in 1:max_iterations
+        changed = false
+        
+        for (i, arg) in enumerate(model_ex.args)
+            if i ∈ trend_equations || !isa(arg, Expr)
+                continue
+            end
+            
+            # For equations LHS = RHS, try to infer degrees
+            if arg.head == :(=)
+                lhs_deg = infer_expression_degree(arg.args[1], degrees, trend_var_set)
+                rhs_deg = infer_expression_degree(arg.args[2], degrees, trend_var_set)
+                
+                # If one side has known degree, propagate to unknowns on other side
+                if lhs_deg !== nothing && rhs_deg === nothing
+                    if propagate_degree!(arg.args[2], lhs_deg, degrees, trend_var_set)
+                        changed = true
+                    end
+                elseif rhs_deg !== nothing && lhs_deg === nothing
+                    if propagate_degree!(arg.args[1], rhs_deg, degrees, trend_var_set)
+                        changed = true
+                    end
+                end
+            end
+        end
+        
+        if !changed
+            break
+        end
+    end
+    
+    # Set remaining unknowns to 0 (stationary)
+    for v in non_trend_vars
+        if degrees[v] == -1.0
+            degrees[v] = 0.0
+        end
+    end
+    
+    return degrees
+end
+
+
+"""
+    infer_expression_degree(expr, degrees::Dict{Symbol, Float64}, trend_vars::Set{Symbol})
+
+Infer the homogeneity degree of an expression given known variable degrees.
+Returns nothing if degree cannot be determined (contains unknowns).
+"""
+function infer_expression_degree(expr, degrees::Dict{Symbol, Float64}, trend_vars::Set{Symbol})
+    if expr isa Number
+        return 0.0
+    elseif expr isa Symbol
+        # Parameter or constant - degree 0
+        return 0.0
+    elseif expr isa Expr
+        if expr.head == :ref
+            var_name = expr.args[1]
+            if var_name isa Symbol
+                deg = get(degrees, var_name, 0.0)
+                return deg == -1.0 ? nothing : deg
+            end
+            return 0.0
+        elseif expr.head == :call
+            op = expr.args[1]
+            
+            if op == :*
+                total = 0.0
+                for i in 2:length(expr.args)
+                    d = infer_expression_degree(expr.args[i], degrees, trend_vars)
+                    if d === nothing
+                        return nothing
+                    end
+                    total += d
+                end
+                return total
+            elseif op == :/
+                if length(expr.args) >= 3
+                    d1 = infer_expression_degree(expr.args[2], degrees, trend_vars)
+                    d2 = infer_expression_degree(expr.args[3], degrees, trend_vars)
+                    if d1 === nothing || d2 === nothing
+                        return nothing
+                    end
+                    return d1 - d2
+                end
+            elseif op == :^
+                base_d = infer_expression_degree(expr.args[2], degrees, trend_vars)
+                if base_d === nothing
+                    return nothing
+                end
+                exp_val = expr.args[3]
+                if exp_val isa Number
+                    return base_d * exp_val
+                else
+                    # Symbolic exponent - assume base degree preserved (for params like α)
+                    return base_d
+                end
+            elseif op in (:+, :-)
+                # For sums, all terms should have same degree
+                if length(expr.args) < 2
+                    return 0.0
+                end
+                first_d = infer_expression_degree(expr.args[2], degrees, trend_vars)
+                if first_d === nothing
+                    return nothing
+                end
+                for i in 3:length(expr.args)
+                    d = infer_expression_degree(expr.args[i], degrees, trend_vars)
+                    if d === nothing
+                        return nothing
+                    end
+                    # In balanced growth, sums can have terms of same degree
+                    # If degrees differ, take the max (dominant term)
+                    if abs(d - first_d) > 0.5
+                        first_d = max(d, first_d)
+                    end
+                end
+                return first_d
+            elseif op in (:exp, :log, :sin, :cos)
+                # Transcendental functions require degree 0 arguments
+                return 0.0
+            end
+        elseif expr.head == :block
+            # Handle block expressions - return degree of last non-LineNumberNode
+            for i in length(expr.args):-1:1
+                if !(expr.args[i] isa LineNumberNode)
+                    return infer_expression_degree(expr.args[i], degrees, trend_vars)
+                end
+            end
+        end
+    end
+    return 0.0
+end
+
+
+"""
+    propagate_degree!(expr, target_degree::Float64, degrees::Dict{Symbol, Float64}, 
+                      trend_vars::Set{Symbol})
+
+Try to assign degrees to unknown variables in expr to achieve target_degree.
+Returns true if any degrees were updated.
+"""
+function propagate_degree!(expr, target_degree::Float64, degrees::Dict{Symbol, Float64},
+                           trend_vars::Set{Symbol})
+    if !(expr isa Expr)
+        return false
+    end
+    
+    changed = false
+    
+    if expr.head == :ref
+        var_name = expr.args[1]
+        if var_name isa Symbol && haskey(degrees, var_name) && degrees[var_name] == -1.0
+            degrees[var_name] = target_degree
+            changed = true
+        end
+    elseif expr.head == :call
+        op = expr.args[1]
+        
+        if op in (:*, :+, :-, :/)
+            # For multiplicative expressions with one unknown, solve for it
+            unknown_count = 0
+            unknown_idx = 0
+            known_sum = 0.0
+            
+            for i in 2:length(expr.args)
+                d = infer_expression_degree(expr.args[i], degrees, trend_vars)
+                if d === nothing
+                    unknown_count += 1
+                    unknown_idx = i
+                else
+                    if op == :*
+                        known_sum += d
+                    elseif op == :/ && i == 3
+                        known_sum -= d
+                    end
+                end
+            end
+            
+            if unknown_count == 1 && op in (:*, :/)
+                # Can solve for the unknown
+                required_deg = target_degree - known_sum
+                if propagate_degree!(expr.args[unknown_idx], required_deg, degrees, trend_vars)
+                    changed = true
+                end
+            else
+                # Recurse into all args with target degree
+                for i in 2:length(expr.args)
+                    if propagate_degree!(expr.args[i], target_degree, degrees, trend_vars)
+                        changed = true
+                    end
+                end
+            end
+        else
+            # For other ops, recurse
+            for i in 2:length(expr.args)
+                if propagate_degree!(expr.args[i], target_degree, degrees, trend_vars)
+                    changed = true
+                end
+            end
+        end
+    elseif expr.head == :block
+        for arg in expr.args
+            if !(arg isa LineNumberNode)
+                if propagate_degree!(arg, target_degree, degrees, trend_vars)
+                    changed = true
+                end
+            end
+        end
+    end
+    
+    return changed
+end
+
+
+"""
+    detrend_variable_ref(var_ref::Expr, degrees::Dict{Symbol, Float64}, 
+                         trend_var::Symbol, growth_param::Symbol)
+
+Transform a variable reference v[t] to its detrended form.
+
+For a variable v with degree d > 0:
+- v[0] → v[0] * T[0]^d (current period)
+- v[-1] → v[-1] * T[-1]^d (lagged)
+- v[1] → v[1] * T[1]^d (lead)
+- v[ss] → v[ss] (steady state - trend cancels in ratios)
+
+Variables with degree 0 are left unchanged.
+"""
+function detrend_variable_ref(var_ref::Expr, degrees::Dict{Symbol, Float64},
+                              trend_var::Symbol, growth_param::Symbol)
+    if var_ref.head != :ref
+        return var_ref
+    end
+    
+    var_name = var_ref.args[1]
+    time_idx = var_ref.args[2]
+    
+    if !(var_name isa Symbol)
+        return var_ref
+    end
+    
+    degree = get(degrees, var_name, 0.0)
+    
+    # Variables with degree 0 or the trend variable itself are unchanged
+    if degree == 0.0 || var_name == trend_var
+        return var_ref
+    end
+    
+    # Check if this is a steady state reference
+    if time_idx isa Symbol && occursin(r"^(ss|stst|steady|steadystate|steady_state)$"i, string(time_idx))
+        # In steady state, detrended variables are just the variable values
+        # because the trend normalizes to 1 on the balanced growth path
+        return var_ref
+    end
+    
+    # Create trend reference with same time index
+    trend_ref = Expr(:ref, trend_var, time_idx)
+    
+    # For degree 1: v[t] → v[t] * T[t]
+    # For other degrees: v[t] → v[t] * T[t]^d
+    if degree == 1.0
+        return Expr(:call, :*, var_ref, trend_ref)
+    else
+        return Expr(:call, :*, var_ref, Expr(:call, :^, trend_ref, degree))
+    end
+end
+
+
+"""
+    detrend_expression(expr, degrees::Dict{Symbol, Float64}, 
+                       trend_var::Symbol, growth_param::Symbol)
+
+Recursively transform an expression by detrending all variable references.
+"""
+function detrend_expression(expr, degrees::Dict{Symbol, Float64},
+                           trend_var::Symbol, growth_param::Symbol)
+    postwalk(x -> begin
+        if x isa Expr && x.head == :ref
+            return detrend_variable_ref(x, degrees, trend_var, growth_param)
+        end
+        x
+    end, expr)
+end
+
+
+"""
+    transform_trend_equation(eq::Expr, trend_var::Symbol, growth_param::Symbol)
+
+Transform the trend equation T[0] = γ * T[-1] into a normalized form.
+
+The original trend equation states that T grows at rate γ.
+For the detrended system, we transform this to: T[0] / T[-1] = γ
+
+This allows the steady state solver to find γ as the ratio.
+Actually, we set T[ss] = 1 as the normalization and remove the trend equation,
+replacing it with: 1 = γ (which is just a constraint on the growth rate steady state).
+
+Or more precisely, we replace the equation with T[0] = T[-1] * γ in terms of 
+the normalized trend where T[ss] = 1.
+"""
+function transform_trend_equation(eq::Expr, trend_var::Symbol, growth_param::Symbol)
+    # The trend equation T[0] = γ * T[-1] in steady state means T = T * γ
+    # This only has a solution if γ = 1 (no growth in steady state)
+    # OR we normalize T[ss] = 1 and accept that the equation holds trivially
+    
+    # Return the equation unchanged - the detrending of other variables
+    # will make them stationary, and T[ss] will be normalized to 1
+    return eq
+end
+
+
+"""
+    apply_automatic_detrending(model_ex::Expr)
+
+Automatically detect trend variables from model equations.
+
+This function:
+1. Detects trend variables from equations like T[0] = γ * T[-1]
+2. Solves homogeneity constraints to determine which variables might be trending
+3. Returns the original equations (unchanged) and balanced growth information
+
+The actual detrending happens during steady state computation where the trend
+is normalized (T[ss] = 1 for deterministic steady state).
+
+Note: This follows the approach in the paper where the model is analyzed for
+balanced growth structure, but the equations are kept in their original form.
+The solver handles the normalization internally.
+"""
+function apply_automatic_detrending(model_ex::Expr)
+    # Step 1: Detect trend variables
+    trend_vars = detect_trend_variables(model_ex)
+    
+    if isempty(trend_vars)
+        # No trends detected - return original equations
+        return model_ex, BalancedGrowthPath()
+    end
+    
+    # Collect all variables from the model
+    all_vars = Set{Symbol}()
+    for arg in model_ex.args
+        if isa(arg, Expr)
+            union!(all_vars, extract_variables_from_expr(arg))
+        end
+    end
+    
+    # Step 2: Solve homogeneity constraints to identify trending variables
+    params = Set{Symbol}()  # Parameters will be identified separately
+    degrees = solve_homogeneity_constraints(model_ex, trend_vars, all_vars, params)
+    
+    # Create balanced growth path info (equations are NOT transformed)
+    growth_params = Set{Symbol}(values(trend_vars))
+    bgp = BalancedGrowthPath(
+        trend_vars,
+        degrees,
+        Expr[],  # No transformed equations - we keep original form
+        growth_params
+    )
+    
+    # Return original equations unchanged - detrending handled in steady state solver
+    return model_ex, bgp
+end
+
+
+"""
+    is_trend_equation(eq::Expr, trend_var::Symbol)
+
+Check if an equation is the trend equation for the given trend variable.
+"""
+function is_trend_equation(eq::Expr, trend_var::Symbol)
+    if eq.head != :(=)
+        return false
+    end
+    
+    lhs = eq.args[1]
+    if lhs isa Expr && lhs.head == :ref && lhs.args[1] == trend_var && lhs.args[2] == 0
+        return true
+    end
+    
+    return false
+end
+
+
+"""
+    analyze_balanced_growth(model_ex::Expr, params::Set{Symbol})
+
+Analyze a model for balanced growth path structure.
+
+Returns a BalancedGrowthPath struct containing:
+- Detected trend variables and their growth rates
+- Homogeneity degrees of all variables
+- Information needed for detrended steady state computation
+"""
+function analyze_balanced_growth(model_ex::Expr, params::Set{Symbol})
+    # Step 1: Detect trend variables
+    trend_vars = detect_trend_variables(model_ex)
+    
+    if isempty(trend_vars)
+        return BalancedGrowthPath()
+    end
+    
+    # Collect all variables from the model
+    all_vars = Set{Symbol}()
+    for arg in model_ex.args
+        if isa(arg, Expr)
+            union!(all_vars, extract_variables_from_expr(arg))
+        end
+    end
+    
+    # Step 2: Solve homogeneity constraints
+    degrees = solve_homogeneity_constraints(model_ex, trend_vars, all_vars, params)
+    
+    # Step 3: Create balanced growth path info
+    growth_params = Set{Symbol}(values(trend_vars))
+    
+    return BalancedGrowthPath(
+        trend_vars,
+        degrees,
+        Expr[],  # Detrended equations computed later if needed
+        growth_params
+    )
+end
+
+
+# ============================================================================
+# End Balanced Growth Path Functions
+# ============================================================================
+
+
 """
 $(SIGNATURES)
 Parses the model equations and assigns them to an object.
@@ -124,6 +796,11 @@ macro model(𝓂,ex...)
 
     model_ex = parse_occasionally_binding_constraints(model_ex::Expr, max_obc_horizon = max_obc_horizon)::Expr
     
+    # Automatic balanced growth path detection and detrending
+    # This step identifies trend variables and transforms the model equations
+    # to work with detrended (stationary) variables
+    model_ex, balanced_growth_info = apply_automatic_detrending(model_ex)
+    
     # obc_shock_bounds = Tuple{Symbol, Bool, Float64}[]
 
     # write down dynamic equations and add auxiliary variables for leads and lags > 1
@@ -841,6 +1518,8 @@ macro model(𝓂,ex...)
                         false, # precompile - to be set by @parameters
 
                         Dict{Symbol, Float64}(), # guess
+                        
+                        $balanced_growth_info, # balanced_growth - detected during parsing
 
                         sort($aux),
                         sort(collect($aux_present)), 
diff --git a/src/structures.jl b/src/structures.jl
index f918b6a37..0855c6e31 100644
--- a/src/structures.jl
+++ b/src/structures.jl
@@ -281,6 +281,35 @@ mutable struct solver_parameters
     backtracking_order::Int
 end
 
+
+"""
+    BalancedGrowthPath
+
+Stores information about the balanced growth path structure of a model.
+
+# Fields
+- `trend_variables::Dict{Symbol, Symbol}`: Maps trend variable names to their growth rate parameters
+  (e.g., `:A => :γ` means A grows at rate γ, i.e., `A[0] = γ * A[-1]`)
+- `variable_degrees::Dict{Symbol, Float64}`: The homogeneity degree of each variable relative to the trend
+  (degree 1 means the variable grows at the same rate as the trend)
+- `detrended_equations::Vector{Expr}`: The model equations after detrending transformation
+- `growth_rate_parameters::Set{Symbol}`: Parameters that represent growth rates
+"""
+struct BalancedGrowthPath
+    trend_variables::Dict{Symbol, Symbol}  # trend_var => growth_rate_param
+    variable_degrees::Dict{Symbol, Float64}  # var => homogeneity degree
+    detrended_equations::Vector{Expr}
+    growth_rate_parameters::Set{Symbol}
+end
+
+BalancedGrowthPath() = BalancedGrowthPath(
+    Dict{Symbol, Symbol}(),
+    Dict{Symbol, Float64}(),
+    Expr[],
+    Set{Symbol}()
+)
+
+
 mutable struct ℳ
     model_name::Any
     # SS_optimizer
@@ -294,6 +323,9 @@ mutable struct ℳ
     precompile::Bool
 
     guess::Dict{Symbol, Float64}
+    
+    # Balanced growth path information
+    balanced_growth::BalancedGrowthPath
 
     # ss
     # dynamic_variables::Vector{Symbol}
diff --git a/test/test_balanced_growth_path.jl b/test/test_balanced_growth_path.jl
new file mode 100644
index 000000000..b198f145a
--- /dev/null
+++ b/test/test_balanced_growth_path.jl
@@ -0,0 +1,154 @@
+using MacroModelling
+using Test
+
+@testset "Balanced Growth Path" begin
+    @testset "Model without balanced growth" begin
+        @model RBC_no_growth begin
+            1  /  c[0] = (β  /  c[1]) * (α * exp(z[1]) * k[0]^(α - 1) + (1 - δ))
+            c[0] + k[0] = (1 - δ) * k[-1] + q[0]
+            q[0] = exp(z[0]) * k[-1]^α
+            z[0] = ρ * z[-1] + std_z * eps_z[x]
+        end
+
+        @parameters RBC_no_growth begin
+            std_z = 0.01
+            ρ = 0.2
+            δ = 0.02
+            α = 0.5
+            β = 0.95
+        end
+
+        @test has_balanced_growth(RBC_no_growth) == false
+        
+        bg_info = get_balanced_growth_path_info(RBC_no_growth)
+        @test bg_info.has_balanced_growth == false
+        @test isempty(bg_info.trend_vars)
+        @test isempty(bg_info.deflators)
+        @test isempty(bg_info.detrended_vars)
+    end
+
+    @testset "Deflator parsing in @model" begin
+        # Test that deflator option is correctly parsed
+        @model RBC_with_deflator deflator = Dict(:c => :A, :k => :A) begin
+            A[0] = γ * A[-1]
+            1/c[0] = β * (1/c[1]) * (α * A[1]^(1-α) * k[0]^(α-1) + 1-δ)
+            c[0] + k[0] = A[0]^(1-α) * k[-1]^α + (1-δ)*k[-1]
+            z[0] = ρ * z[-1] + σ * eps[x]
+        end
+
+        @parameters RBC_with_deflator trend_var = Dict(:A => :γ) begin
+            γ = 1.02
+            α = 0.33
+            β = 0.99
+            δ = 0.025
+            ρ = 0.9
+            σ = 0.01
+        end
+
+        @test has_balanced_growth(RBC_with_deflator) == true
+        
+        bg_info = get_balanced_growth_path_info(RBC_with_deflator)
+        @test bg_info.has_balanced_growth == true
+        @test :A ∈ keys(bg_info.trend_vars)
+        @test bg_info.trend_vars[:A] == :γ
+        @test :c ∈ keys(bg_info.deflators)
+        @test :k ∈ keys(bg_info.deflators)
+        @test bg_info.deflators[:c] == :A
+        @test bg_info.deflators[:k] == :A
+        @test :c ∈ bg_info.detrended_vars
+        @test :k ∈ bg_info.detrended_vars
+    end
+
+    @testset "Trend_var parsing in @parameters" begin
+        # Test that trend_var is correctly stored when specified in @parameters
+        # Use the same model structure as previous test to avoid analytical solve edge cases
+        @model TrendVarTest deflator = Dict(:y => :A) begin
+            A[0] = γ * A[-1]
+            1/c[0] = β * (1/c[1]) * (α * y[1]/k[0] + 1-δ)
+            y[0] = k[-1]^α
+            c[0] + k[0] = y[0] + (1-δ)*k[-1]
+        end
+
+        @parameters TrendVarTest trend_var = Dict(:A => :γ) begin
+            γ = 1.015
+            α = 0.3
+            β = 0.99
+            δ = 0.1
+        end
+
+        bg_info = get_balanced_growth_path_info(TrendVarTest)
+        @test :A ∈ keys(bg_info.trend_vars)
+        @test bg_info.trend_vars[:A] == :γ
+    end
+
+    @testset "Empty deflator still works" begin
+        # Test that model works when deflator is explicitly empty
+        @model RBC_empty_deflator deflator = Dict{Symbol,Symbol}() begin
+            1  /  c[0] = (β  /  c[1]) * (α * exp(z[1]) * k[0]^(α - 1) + (1 - δ))
+            c[0] + k[0] = (1 - δ) * k[-1] + q[0]
+            q[0] = exp(z[0]) * k[-1]^α
+            z[0] = ρ * z[-1] + std_z * eps_z[x]
+        end
+
+        @parameters RBC_empty_deflator begin
+            std_z = 0.01
+            ρ = 0.2
+            δ = 0.02
+            α = 0.5
+            β = 0.95
+        end
+
+        @test has_balanced_growth(RBC_empty_deflator) == false
+    end
+
+    @testset "Automatic trend detection - detect_trend_variables function" begin
+        # Test the detect_trend_variables function directly
+        # Pattern: X[0] = γ * X[-1]
+        model_ex1 = :(begin
+            A[0] = γ * A[-1]
+            y[0] = A[0] * k[-1]^α
+        end)
+        
+        trends1 = MacroModelling.detect_trend_variables(model_ex1)
+        @test :A ∈ keys(trends1)
+        @test trends1[:A] == :γ
+        
+        # Pattern: X[0] = X[-1] * γ (reversed order)
+        model_ex2 = :(begin
+            B[0] = B[-1] * μ
+            y[0] = B[0] * n[0]
+        end)
+        
+        trends2 = MacroModelling.detect_trend_variables(model_ex2)
+        @test :B ∈ keys(trends2)
+        @test trends2[:B] == :μ
+        
+        # No trend variables in standard RBC
+        model_ex3 = :(begin
+            1/c[0] = β * (1/c[1]) * (α * k[0]^(α-1) + 1-δ)
+            c[0] + k[0] = k[-1]^α + (1-δ)*k[-1]
+            z[0] = ρ * z[-1] + σ * eps[x]
+        end)
+        
+        trends3 = MacroModelling.detect_trend_variables(model_ex3)
+        @test isempty(trends3)
+    end
+
+    @testset "Automatic deflator detection - detect_deflators_from_equations function" begin
+        # Test the detect_deflators_from_equations function directly
+        model_ex = :(begin
+            A[0] = γ * A[-1]
+            y[0] = A[0] * k[-1]^α
+            c[0] + k[0] = y[0] + (1-δ)*k[-1]
+        end)
+        
+        trend_vars = Dict{Symbol, Union{Symbol, Expr, Number}}(:A => :γ)
+        deflators = MacroModelling.detect_deflators_from_equations(model_ex, trend_vars)
+        
+        # Variables in the same equations as A should be assigned A as deflator
+        @test :y ∈ keys(deflators)
+        @test :k ∈ keys(deflators)
+        @test deflators[:y] == :A
+        @test deflators[:k] == :A
+    end
+end