Question: How to analytically compute omegabar steady-state value in MacroModelling

[Wastonvv]

Hi, thanks for this excellent package — it’s been incredibly helpful for building and solving DSGE models efficiently!

I’m currently implementing the BGG (1999) model using MacroModelling.jl, and I encountered some difficulty in defining the steady-state value of omegabar.

In this model, omegabar is implicitly defined through the following nonlinear condition:

function comp_omegabar(x, sigma, mu, L)
    F = cdf(Normal(), log(x)/sigma + sigma/2)
    F_p = 1/(sigma*x) * pdf(Normal(), log(x)/sigma + sigma/2)
    G = cdf(Normal(), log(x)/sigma - sigma/2)
    G_p = 1/sigma * pdf(Normal(), log(x)/sigma + sigma/2)
    Gamma = G + x * (1 - F)
    Gamma_p = 1 - F
    err = (Gamma_p - mu*G_p)/Gamma_p * (1 - Gamma) - (Gamma - mu*G)/(L - 1)
    return err
end

And in the @parameters block I tried to define it as:

omegabar[ss] -> fzero(x -> comp_omegabar(x, sigm, mu, L[ss]), 0, 2)

However, it seems that MacroModelling cannot directly handle this definition — it either fails to evaluate at steady state or gets stuck before computing the full model solution.

here is the complete code

@model BGG_1999 begin

    1/C[0] - lambd[0]=0

    zeta/M[0] - lambd[0] + bet *lambd[0] /pi1[1] = 0

    bet * lambd[1] *R[0] = lambd[0]

    lambd[0] *w[0] = xi /(1-H[0] )

    F[0] = normcdf(log(omegabar[0] )/sigm + sigm/2)
    G[0] = normcdf(log(omegabar[0] )/sigm - sigm/2)
    Gamm[0] = G[0] + omegabar[0] * (1 - F[0] )
    S[0] = Rk[1] /R[0] 

    L[0] = Q[0] * K[0] /N[0]                   
    
    (Gammap[0] - mu * Gp[0] ) / Gammap[0] *(1-Gamm[0] ) + (Gamm[0] - mu *G[0] )=1/S[0]
    (Gamm[0] -mu*G[0] )*S[0] = 1 - 1 / L[0]
    V[0] =(1 -  Gamm[0]) * S[0] *R[0] *Q[0] *K[0]        
    Ce[0] = (1-gamma_e)*V[0]
    N[0] = gamma_e * V[-1] + He * w_e[-1] # 
    Ye[0] =A[0] *K[-1] ^alph *(H[0] ^ (Omega) * He ^(1-Omega))^(1-alph) # 
    w[0] * H[0] =(1-alph) * Omega * pm[0] * Ye[0]
    w_e[0] * He =(1-alph) * (1-Omega) * pm[0] *Ye[0]
    Rk[0] = alph * pm[0] *Ye[0] / Q[-1] /K[-1] + Q[0] *(1-delt)/Q[-1] 
    pistar[0] = epsilon/(epsilon-1)* x1[0] /x2[0]
    x1[0] = lambd[0] * pm[0] * Y[0] + thet*bet*x1[1] *(pi1[1])^epsilon
    x2[0] =lambd[0] * Y[0] + thet*bet * x2[1] *(pi1[1])^(epsilon-1)
    thet * pi1[0] ^(epsilon-1)+(1-thet)* pistar[0] ^(1-epsilon) = 1
    K[0] =(1-delt) * K[-1] +I[0] - chi/2*(I[0] /K[-1] -delt)^2*K[-1] 
    Q[0] =1/(1-chi * (I[0] /K[-1] - delt))  
    Ye[0] = D[0] *Y[0]  
    D[0] = thet * pi1[0] ^epsilon*D[-1] +(1-thet) *pistar[0] ^(-epsilon)
    Y[0] =C[0] +Ce[0] +I[0] +chi/2 *(I[0] /K[-1] -delt)^2*K[-1] +mu *G[0] *Rk[0] *Q[-1] *K[-1]  
    pi1[1] * R[0] = In[0] 
    log( In[0] / Ins ) = rho * log( In[-1] / Ins )+ (1-rho)* phipi* log(pi1[0] /pi1s ) + epsr[x]
    log(A[0] )= rhoa * log(A[-1] ) + epsa[x]
    # He[0] = 1
    Gammap[0] = 1 - F[0] 
    Gp[0] = normpdf(log(omegabar[0] )/ sigm + sigm/2) / sigm 
    
end

function comp_omegabar(x,sigma , mu ,L)
    F=cdf(Distributions.Normal() ,(log(x))/sigma+sigma/2 )
    F_p= 1/(sigma*x)*pdf(Distributions.Normal(),(log(x))/sigma+sigma/2);
    G=cdf(Distributions.Normal(),(log(x))/sigma-sigma/2);
    G_p=1/sigma*pdf(Distributions.Normal(),log(x)/sigma+sigma/2);
    Gamma=G+x*(1-F);
    Gamma_p=1-F;
    err=(Gamma_p-mu*G_p)/Gamma_p*(1-Gamma)-(Gamma-mu*G)/(L-1);
    return err
end


@parameters BGG_1999 begin
    sigm=0.48; 
    mu=0.12; 
    bet=0.99; 
    zeta=3.33; 
    alph=0.35;
    chi=1-0.64/(1-alph); 
    delt=0.025; 
    thet=0.8;
    rhoa=0.9;
    phipi=1.5;  
    rho=0.9;
    epsilon=11;
    Omega=0.8;
    xi = 1.4
    He = 1 
    Ins = 1/bet
    pi1s = 1
    gamma_e | L[ss] = 3/2
    omegabar[ss] ->  fzero( x -> comp_omegabar(x, sigm , mu, L[ss]) ,  0, 2) 
end
get_solution(BGG_1999)(:,:omegabar)

ERROR: AssertionError: Could not find stable first order solution.

[Wastonvv]

This setup worked well in Julia 1.9 with MacroModelling.jl v0.1.39,
but it no longer works in Julia 1.11 with MacroModelling.jl v0.1.41

[thorek1]

Hi,

From what I can tell without running the code is that you shouldn't need to define omegbar[ss]. The model equations in steady state should be enough to uniquely pin it down.

You cannot use your own functions inside the model and parameter macros. You shouldn't need to either and in your case you definitely don't need to.

Try it without the omegabar[ss] line in the steady state block.

Let me know if it works.

[Wastonvv]

After further testing, I found something quite puzzling regarding the steady-state computation behavior across different Julia and MacroModelling versions.

When I comment out the line

# omegabar[ss] -> fzero(x -> comp_omegabar(x, sigm, mu, L[ss]), 0, 2)

the model still computes successfully, though quite slowly, under:

Julia 1.10

MacroModelling v0.1.39

Machine: MacBook Pro (M1 Pro chip)

The log output looks like this:

Remove redundant variables in non stochastic steady state problem:      0.894 seconds
Set up non stochastic steady state problem:                             2.676 seconds
Take symbolic derivatives up to first order:                            11.249 seconds
Find non stochastic steady state:                                       114.252 seconds
Model:        BGG_1999
Variables
 Total:       32
  Auxiliary:  0
 States:      7
  Auxiliary:  0
 Jumpers:     5
  Auxiliary:  0
Shocks:       2
Parameters:   18
Calibration
equations:    1

However, using the exact same code under:

Julia 1.11

MacroModelling v0.1.41

Same machine (MacBook Pro, M1 Pro chip)

the solver fails to find a steady state, and gives the following warning:

Remove redundant variables in non-stochastic steady state problem:      1.031 seconds
Set up non-stochastic steady state problem:                             11.429 seconds
Find non-stochastic steady state:                                       122.02 seconds
┌ Warning: Could not find non-stochastic steady state. Consider setting bounds on variables or calibrated parameters in the `@parameters` section (e.g. `k > 10`).
└ @ Main ~/.julia/packages/MacroModelling/WLYku/src/macros.jl:1558
Take symbolic derivatives up to first order:                            3.797 seconds
Model:        BGG_1999
Variables
 Total:       32
  Auxiliary:  0
 States:      7
  Auxiliary:  0
 Jumpers:     5
  Auxiliary:  0
Shocks:       2
Parameters:   18
Calibration
equations:    1

Even more interestingly, when I tested the same code on a Windows machine running
Julia 1.9 + MacroModelling v0.1.39, the steady state was found very quickly (first run ≈7.8 seconds).

So it seems that:

Removing the omegabar[ss] definition is not crucial for solving,

but the performance and convergence behavior varies drastically across Julia and package versions, and possibly by platform (macOS ARM vs Windows).

Could this be related to recent changes in how MacroModelling handles nonlinear steady-state solvers or symbolic differentiation in Julia ≥1.11?

[thorek1]

thanks for testing this out.

I just checked myself and what happens is that he finds the steady state using simulated annealing over the solver parameters (see here for details - item 5 third bullet on the slide).

this approach is stochastic (replicate with setting a seed) and does not guarantee finding a steady state but it significantly improves the chances.

I noticed that the package sometimes throws a warning after the @parameter macro even though if you then call SS(BGG_1999) he does return a valid steady state.

I will look into that, and make sure the warning is only thrown if he really doesn't find the steady state.

as for the first order solution the model appears to have other issues (check the steady state, equations, and parameter values if they make sense)

[Wastonvv]

Thanks a lot for the clarification. I’ll look into the simulated annealing approach .

Really appreciate your detailed explanation and all your work on improving MacroModelling! 🙏