ERROR: LoadError: UndefVarError: x not defined

[gpetrini]

Hi,

I am trying to replicate the model of this paper here which was originally written in matlab (attached as txt, but it is .m). I am not proficient in neither matlab, julia, or even DSGE. As a test, I replicate the RBC model (under RBC folder) and achieve the same results as the matlab code. However, I got the following error when running the baseline model: ERROR: LoadError: UndefVarError: μᶜ not defined

I tried to run a debugger and expand the macros, but no success in finding the problem. Do you have suggestions to debug it?

original matlab code

PS1: the evaluation is considerably slow in comparisson with other codes that I tried PS2: the original paper uses AIM algorith to solve the model PS3: Thanks for the package :)

Model

Here is the complete code of the model:

@model Baseline begin
    ## Resource constraint
    ynet[0] = (c_yc/(ynetmcons)) * c[0] + (sk_yc/(ynetmcons)) * sk[0] + (sc_yc/(ynetmcons)) * sc[0] + actualhc_yc/(ynetmcons) * hc[0] + (ac_zc/(1-ac_zc)*actualhc_yc/(ynetmcons)) * zc[0] - (inv(zc_ac-1)*actualhc_yc/(ynetmcons)) * ac[0] + (actualhk_yc/(ynetmcons)) * hk[0] + (ak_zk/(1-ak_zk)*actualhk_yc/(ynetmcons)) * zk[0] - (inv(zk_ak-1)*actualhk_yc/(ynetmcons)) * ak[0]
    # ## Euler equation
    c[0] = c[1] - r[0]
    # ## Aggregate production function FIXME 0 or ss?
    yc[0] = (1/(1-bb)) * chiz[0] + α * k[-1] - α * l[0] + ly[0] + ((μᶜ[0] - 1)/(1 - bb)) * Nᶜ[0] - ((bb+μᶜ[0]*log(Nᶜ[0]))/(1-bb)) * μᶜ[0] + α * u[0] + (bb/(1-bb)*(θ-1)) * ac[0]
    # ## demand of capital
    yc[0] = kl[0] - l[0] + ly[0] + μᶜ[0] + d[0] * (1/(1+δ/d_pk)) + pk[0] * ((δ/(d_pk+δ))) + u[0] * edu*(1/(d_pk/δ+1))
    # ## capacity choice
    yc[0] = (1 + edp) * u[0] + μᶜ[0] + kl[0] - l[0] + ly[0] + pk[0]
    # ## labor demand in final output
    yc[0] = μᶜ[0] + ly[0] + w[0]
    # ## production of new investment goods FIXME 0 or ss?
    yk[0] = chi[0] * (1 / ( 1 - bb)) + α * kl[0] - l[0] * (α-1/(1-lc_l)) + u[0] * α - ly[0] * (lc_l/(1-lc_l)) + Nᵏ[0] * ((μᵏ[0]-1)/(1-bb)) - ((bb+μᵏ[0]*log(Nᵏ[0]))/(1-bb)) * μᵏ[0] + pk[0] * (bb/(1-bb)) + ak[0] * (bb/(1-bb)*(θ-1))
    # ## Real wage
    w[0] = ζ * l[0] + c[0] + μʷ[0]
    # ## Profits embodied
    prk[0] = yk[0] + pk[0] - μᵏ[0]
    # ## Profits disembodied
    prc[0] = yc[0] - μᶜ[0]
    # ## Value of an adopted innovation for embodied
    vk[0] = ak[0] * ((1-prk_vk)) + prk[0] * prk_vk + vk[1] * (1 - prk_vk)
    # ## Value of an adopted innovation for disembodied
    vc[0] = ac[0] * (1-prc_vc) + prc[0] * (prc_vc) + vc[1] * (1 - prc_vc) - ac[1] * ((1-prc_vc)) - r[0] * ((1-prc_vc))
    # ## Capital accumulation
    k[0] = - u[0] * (edu*δ/(1+gk)) + k[-1] * jcof + yk[0] * (1 - jcof)
    # ## Law of motion for embodied productivity
    zk[0] = zk[-1] + sk[-1]  * (ρ*(gzk+ok)/(1+gzk)) - sf[1] * (ρ*(gzk+ok)/(1+gzk)) + chik[1] * ((gzk+ok)/(1+gzk))
    # ## Law of motion for disembodied productivity
    zc[0] = zc[-1] + sc[1] * (ρ*(gzc+oc)/(1+gzc)) -  sf[1] * (ρ*(gzc+oc)/(1+gzc))
    # ## Free entry for embodied
    sk[0] * (1 - ρ) = zk[0] - sf[0] * (ρ) + jk[1] - zk[1] - r[0]
    # ## Free entry for disembodied
    sc[0] * (1 - ρ) = zc[0] - sf[0] * (ρ) + jc[1] - zc[1] - r[0]
    # ## Bellman for not adopted disemb innovation
    - jc[0] = hc[0] * (hc_jc+ϕᶜ*elc*λᶜ/R*rz*(1-zc_ac*vc_jc)) + r[0] * ((1+hc_jc)) - zc[0] * (ϕᶜ*rz*((1-λᶜ)+λᶜ*zc_ac*vc_jc)/R) + ac[1] * (ϕᶜ*λᶜ*rz*zc_ac*vc_jc/R) - vc[1] * (ϕᶜ*λᶜ*rz*zc_ac*vc_jc/R) + sf[0] * (ϕᶜ*elc*λᶜ*rz/R*(zc_ac*vc_jc-1)) + zc[1] * (ϕᶜ*rz*(1-λᶜ)/R) - jc[1] * (ϕᶜ*rz*(1-λᶜ)/R)
    # ## law of motion for adopted disembodied innvo
    ac[0] = ac[-1] * (ϕᶜ*(1-λᶜ)/(1+gzc)) + hc[-1] * (elc*λᶜ*((ϕᶜ/(1+gzc))*zc_ac-ϕᶜ/(1+gzc))) - sf[-1] * (elc*λᶜ*((ϕᶜ/(1+gzc))*zc_ac-ϕᶜ/(1+gzc))) + zc[-1] * ((1-ϕᶜ*(1-λᶜ)/(1+gzc)))
    # ## optimal investment in adoption of disemb innov
    zc[0] = sf[0] * (1+ℓᶜ) + r[0] - hc[0] * ℓᶜ - vc[1] * (1/(1-jc_vc*ac_zc)) + ac[1] * (1/(1-jc_vc*ac_zc)) + jc[1] * (1/(vc_jc*zc_ac-1)) - zc[1] * (1/(vc_jc*zc_ac-1))
    # ## Bellman for not adopted emb innovation
    - jk[0] = hk[0] * ((hk_jk+(1-ok)*elk*λᵏ/R*ra*(1-zk_ak*vk_jk))) + r[0] * ((1+hk_jk)) - zk[0] * (ϕᵏ*ra*((1-λᵏ)+λᵏ*zk_ak*vk_jk)/R) + ak[1] * (ϕᵏ*λᵏ*ra*zk_ak*vk_jk/R) - vk[1] * (ϕᵏ*λᵏ*ra*zk_ak*vk_jk/R) + sf[0] * (ϕᵏ*elk*λᵏ*ra/R*(zk_ak*vk_jk-1)) + zk[1] * (ϕᵏ*ra*(1-λᵏ)/R) - jk[1] * (ϕᵏ*ra*(1-λᵏ)/R)
    # ## law of motion for adopted embodied innvo
    ak[0] = ak[-1] * (ϕᵏ*(1-λᵏ)/(1+gzk)) + hk[-1] * (elk*λᵏ*((ϕᵏ/(1+gzk))*zk_ak-ϕᵏ/(1+gzk))) - sf[-1] * (elk*λᵏ*((ϕᵏ/(1+gzk))*zk_ak-ϕᵏ/(1+gzk))) + zk[-1] * ((1-ϕᵏ*(1-λᵏ)/(1+gzk)))
    # ## optimal investment in adoption of emb innov
    zk[0] = sf[0] * ((1+ℓᵏ)) + r[0] - hk[0] * (ℓᵏ) - vk[1] * (1/(1-jk_vk*ak_zk)) + ak[1] * (1/(1-jk_vk*ak_zk)) + jk[1] * (1/(vk_jk*zk_ak-1)) - zk[1] * (1/(vk_jk*zk_ak-1))
    # ## Arbitrage
    pk[0] = -r[0] + d[1] * ((R-1-gpk)/R) + pk[1] * ((1+gpk)/R)
    # ## entry into final goods sector
    μᶜ[0] = - yc[0] * mucof + sf[0] * mucof + Nᶜ[0] * mucof
    # ## m
    μᶜ[0] = Nᶜ[0]  * ημᶜ
    # ## entry into capital goods sector
    μᵏ[0] = - mukcof * yk[0] - pk[0] * mukcof + sf[0] * mukcof + Nᵏ[0] * mukcof
    # ## mk
    μᵏ[0] = Nᵏ[0] * ημᵏ
    # ## equivalence between klzero and jlag
    kl[0] = k[-1]
    # ## Definition of output net of total overhead costs
    ynet[0] = yc[0] * (1/(1-oc_yc)) - Nᶜ[0] * (occ_yc/(1-oc_yc)) - Nᵏ[0] * (ock_yc/(1-oc_yc)) - sf[0] * (oc_yc/(1-oc_yc))
    # ## definition of scaling factor
    sf[0] = kl[0] + ak[0] * (bb*(1-θ)) - ac[0] * (bb*(1-θ))
    # ## definition of ynetm
    ynetm[0] = ynet[0] * (1/(1-mc_yc*inv(ynet_yc)-mk_yc*inv(ynet_yc))) - yc[0] * (mc_yc/ynetmcons) + μᶜ[0] * (mc_yc/ynetmcons) - pk[0] * (mk_yc/ynetmcons) - yk[0] * (mk_yc/ynetmcons) + μᵏ[0] * (mk_yc/ynetmcons)
    # ## Definition of total value added
    yT[0] = ynetm[0] * (ynetmcons/(ynetmcons+pkyk_yc)) + pk[0] * (pkyk_yc/(ynetmcons+pkyk_yc)) + yk[0] * (pkyk_yc/(ynetmcons+pkyk_yc))
    # ## labor demand in capital goods production
    yk[0] = - pk[0] + μᵏ[0] + w[0] + l[0] * (1/(1-lc_l)) - ly[0] * (lc_l/(1-lc_l))
    # ## embodied productivity shock process
    chi[0] = ρᵡ * chi[-1] + σᵡ* eps_chi[x]
    # ## Labor augmenting technology shock process
    chiz[0] = ρᶻᵪ * chiz[-1] + σᶻᵪ * eps_chi_z[x]
    # ## Wage markup shock process
    μʷ[0] = μʷ[-1] * ρᵐʷ + σᵐʷ * eps_muw[x]
    # ## Wage markup shock process
    chik[0] = ρᵏᵪ * chik[-1] + σᵏᵪ * eps_chi_k[x]

end

@parameters Baseline simplify = false guess = Dict(:μʷ => 1.2, :μᶜ => 1.1, :μᵏ => 1.2, :Nᶜ => 1, :Nᵏ => 1, :u => 0.8)  begin

    β = 0.95 # discount factor (bet)
    δ = 0.08 # depreciation rate
    ζ = 1 # labor supply curvature (fi)
    α = 1/3 # k share (al)
    g_y = 0.2 * 0.7 # ss g/y ratio
    ν = 0.6 # FIXME Comment
    θ = 1 / ν # elasticity of substitution intermediate good sector (th)
    ρ = 0.9 # parameter embodied technology
    η = 0.0
    # μʷ[ss] = 1.2 # ss wage markup (muw)
    # μᶜ[ss] = 1.1 # ss wage markup (muc)
    # Nᶜ[ss] = 1 # (nc)
    # Nᵏ[ss] = 1 # (nk)
    dμᶜ = μᶜ[ss] # (dmuc)
    ημᶜ = (dμᶜ / μᶜ[ss]) * Nᶜ[ss]
    boc = (μᶜ[ss] - 1)/μᶜ[ss]
    # μᵏ[ss] = 1.2
    ημᵏ = ημᶜ
    λᵏ = 0.1 # (lamk)
    λᶜ = 0.1 # (lamc)
    elk = 0.9
    elc = 0.9
    ℓᵏ = elk - 1 # elasticity of λ wrt H
    ℓᶜ = elc - 1 # elasticity of λ wrt H
    o = 0.03
    oz = 0.03
    oc = 0.03
    ok = 0.03
    ϕᶜ = 1 - oc # (phic)
    ϕᵏ = 1 - ok # (phic)
    bb = 0.5 # intermediate share in final output

    ## Nonstochastic steady state
    gpk    = -0.026
    gy     =  0.024
    gk     =  gy - gpk
    gzc    = (gy-α*gk)/bb*(1-bb)/(θ-1)
    gzk    = (gpk-gzc*bb*(θ-1))/(bb*(1-θ))
    gtfp   =  gy-α*gk+gzk*(α*bb*(θ-1))/(1-α*(1-bb))
    measbls = (0.014-gy+α*gk)/(gzk*(α*bb*(θ-1))/(1-α*(1-bb)))
    gv     = gy
    gvz    = gy
    R      = (1+gy)/β
    d_pk   = R-(1+gpk)                   # definition of R
    yc_pkkc = μᶜ[ss]/(α*(1-bb))*(d_pk+δ) # foc for k
    yk_kk   = μᵏ[ss]/(α*(1-bb))*(d_pk+δ) # new capital to capital in capital production sector
    yk_k   = (gk+δ)/(1+gk)             # new capital to capital ratio
    kk_k   = yk_k/yk_kk                  # share of capital in capital production.
    kc_k   = 1-kk_k
    kk_kc  = kk_k/kc_k
    lk_lc  = kk_kc
    lk_l   = lk_lc/(lk_lc+1)
    lc_l   = 1-lk_l
    pkyk_yc= kk_kc*μᵏ[ss]/μᶜ[ss]
    mk_yc  = bb*1/θ*pkyk_yc/μᵏ[ss]
    mc_yc  = bb*1/θ/μᶜ[ss]
    pkk_yc = inv(yc_pkkc)/kc_k
    pik_yc = pkk_yc*μᶜ[ss]/μᵏ[ss]              # value of total capital stock removing fluctuations in relative price of capital due to markup variations
    prk_yc   = pkyk_yc*(1-1/θ)*bb/μᵏ[ss]
    prc_yc  = (1-1/θ)*bb/μᶜ[ss]
    prk_vk   = 1-(1+gv)*ϕᵏ/((1+gzk)*R) # bellman for va
    prc_vc = 1-(1+gvz)*ϕᶜ/((1+gzc)*R) # bellman for vz
    yc_vk    = prk_vk*inv(prk_yc)
    yc_vc   = prc_vc*inv(prc_yc)
    zk_ak   = ((gzk+ok)/(λᵏ*ϕᵏ)+1)
    zc_ac   = ((gzc+oc)/(λᶜ*ϕᶜ)+1)
    ac_zc   = inv(zc_ac)
    ak_zk   = inv(zk_ak)
    ra     = (1+gy)/(1+gzk)
    rz     = (1+gy)/(1+gzc)

    jk_yc    = (inv(1-elk*ϕᵏ*λᵏ*ra/R-(1-λᵏ)*ϕᵏ*ra/R))*(1-elk)*ϕᵏ*λᵏ*ra*zk_ak/R*inv(yc_vk) # zk * jk /yc bellman for not adopted innov
    jc_yc   = (inv(1/ϕᶜ-elc*λᶜ*rz/R-(1-λᶜ)*rz/R))*(1-elc)*λᶜ*rz*zc_ac/R*inv(yc_vc) # zc*jc/yc bellman for not adopted innov

    hk_yc    = ϕᵏ*elk*λᵏ*ra/R*(inv(yc_vk)*zk_ak-jk_yc) # zk *hk/yc
    hc_yc    = ϕᶜ*elc*λᶜ*rz/R*(inv(yc_vc)*zc_ac-jc_yc) # zc *hc/yc

    sk_yc    = jk_yc*(gzk+o)*(1+gv)*inv((1+gzk)*R) # from free entry cond't
    sc_yc   = jc_yc*(gzc+oz)*(1+gvz)*inv((1+gzc)*R)

    hc_jc  = hc_yc/jc_yc
    hk_jk  = hk_yc/jk_yc

    vc_jc  = inv(yc_vc)/jc_yc
    vk_jk  = inv(yc_vk)/jk_yc

    jc_vc=inv(vc_jc)
    jk_vk=inv(vk_jk)

    bock   = boc*pkyk_yc*(μᵏ[ss]-1)*μᶜ[ss]/(μᵏ[ss]*(μᶜ[ss]-1))
    occ_yc=boc*pik_yc
    ock_yc=bock*pik_yc
    oc_yc=occ_yc+ock_yc

    c_yc = 1-oc_yc-g_y-mc_yc-mk_yc-sk_yc-sc_yc-((ϕᶜ/(1+gzc))^2-inv(zc_ac))*hc_yc-((ϕᵏ/(1+gzk))^2-inv(zk_ak))*hk_yc
    pi_yc=(μᶜ[ss]-1)/μᶜ[ss]-oc_yc
    # u[ss]=.8
    edu=α*(1-bb)*yc_pkkc/(μᶜ[ss]*δ) # from foc wrt utilization, edu = elasticity of depreciation with respect to capacity
    edup=0 # partial of edu wrt u
    edp=1/3 # (edu)-1+edup/(edu*u) # elasticity of del' (i.e. elasticity of delta prima)
    actualhk_yc=hk_yc*(1-ak_zk) # total expenses in adoption of capital specific innovations
    actualhc_yc=hc_yc*(1-ac_zc) # total expenses in adoption of consumption specific innovations
    inv_Y=pkyk_yc/(pkyk_yc+1-mc_yc-mk_yc-occ_yc-ock_yc) # investment output ratio
    Y_yc=pkyk_yc/inv_Y
    ynet_yc=1-oc_yc

    # NOTE Coefficients for the log-linearization
    qcof = (1-δ)*(1+gpk)/R
    jcof = (1-δ)/(1+gk)
    vcof = (1+gy)/((1+gzk)*R)
    vzcof= (1+gy)/((1+gzc)*R)
    mucof= μᶜ[ss]-1
    mukcof=μᵏ[ss]-1

    ycof1=ynet_yc*inv(ynet_yc-mc_yc-mk_yc+pkyk_yc-actualhk_yc-actualhc_yc-sk_yc-sc_yc)
    ycof2=mc_yc*inv(ynet_yc-mc_yc-mk_yc+pkyk_yc-actualhk_yc-actualhc_yc-sk_yc-sc_yc)
    ycof3=mk_yc*inv(ynet_yc-mc_yc-mk_yc+pkyk_yc-actualhk_yc-actualhc_yc-sk_yc-sc_yc)
    ycof4=pkyk_yc*inv(ynet_yc-mc_yc-mk_yc+pkyk_yc-actualhk_yc-actualhc_yc-sk_yc-sc_yc)
    ycof=inv(ynet_yc-mc_yc-mk_yc+pkyk_yc-actualhk_yc-actualhc_yc-sk_yc-sc_yc)

    ynetmcons=1-oc_yc-mc_yc-mk_yc # fraction of ynetm in y

    # NOTE Shock to Embodied Technology
    ρᵡ = (0.7)^4   # autoregressive component
    σᵡ = 0.01    # standard deviation

    # Disembodied Technology Shock
    ρᶻᵪ = (0.7)^4
    σᶻᵪ = 0.01

    # Wage markup shock
    ρᵐʷ   = 0.60 # (rhow)
    σᵐʷ = 0.01

    # Chik shock
    ρᵏᵪ = 0.8
    σᵏᵪ = 0.01

end

Log message

Take symbolic derivatives up to first order:                1598.541 seconds
Find non stochastic steady state:                   ERROR: LoadError: UndefVarError: `μᶜ` not defined
Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/MacroModelling/oKvYN/src/MacroModelling.jl:3501 [inlined]
  [2] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:163 [inlined]
  [3] macro expansion
    @ ./none:0 [inlined]
  [4] generated_callfunc
    @ ./none:0 [inlined]
  [5] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Vector{…}, ::MacroModelling.ℳ, ::Bool, ::Bool, ::Vector{…})
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:150
  [6] macro expansion
    @ ~/.julia/packages/MacroModelling/oKvYN/src/macros.jl:1518 [inlined]
  [7] top-level scope
    @ ~/Documents/KS_vs_DSGE/code/DSGE/Comin_Gertler_2006/Baseline.jl:1513
  [8] include(fname::String)
    @ Base.MainInclude ./client.jl:489
  [9] top-level scope
    @ REPL[1]:1
 [10] top-level scope
    @ none:1
in expression starting at /home/gpetrini/Documents/KS_vs_DSGE/code/DSGE/Comin_Gertler_2006/Baseline.jl:102

[thorek1]

Hi,

thank you for the posting here.

I had a brief look and there are multiple equations of the sorts: dμᶜ = μᶜ[ss]

This does not work because you make a variable dependent on a steady state value which then influences other or even the same steady state value (via ημᵏ in this case). This system of equations is underdetermined. In order to fix that you can endogenously determine a parameter value based on a calibration relationship. Have a look at the @parameters docstring for examples. This is certainly more involved compared to just fixing these parameters to steady state values as you can with other software. Doing it the suggested way means you have more relationships between the parameters. I'll think about writing an interface for the other way.

Another way to get you going is to explicitly hard code the parameter values without the references to steady state values (e.g. you know μᶜ[ss] = 1/3 so set dμᶜ = 1/3).

@PS1: let me know about the runtime once you fixed the steady state issues @PS2: it appears most time is spent taking derivatives not solving the model. this is hopefully resolved with the steady state issue @PS3: you are welcome. let me know if you run into further roadblocks and I will see if I can find some time to write better error messages for the model setup

[gpetrini]

Thank you a lot for your support!

I followed the second approach, but I am planning to use your suggestion based on the calibration relationship.

The code is now running. I cannot replicate the original paper so far, but there is certainly some error in my script.

One issue that I am unable to solve is the performance of the code. It takes 1400 s (at best) to finish the symbolic derivatives (up to 1st order). Do you have any advice on how to speed up the computation? (or to spot bottlenecks?).

The model is not particularly huge (and there are larger models in the examples folder), so I imagine that I am introducing some slowness in it.

Thank you once again!

[thorek1]

Glad to hear it works.

1400s sounds way too long indeed. I can have a look. Could you post the working model code? Which Julia and MacroModelling version are you using?

[gpetrini]

Sure!

Thank you in advance.

Julia and MacroModelling version

(@v1.10) pkg> status
Status `~/.julia/environments/v1.10/Project.toml`
  [00ebfdb7] CSTParser v3.4.3
  [a93c6f00] DataFrames v1.6.1
  [31a5f54b] Debugger v0.7.9
  [5203de40] Dynare v0.9.16
  [5903a43b] Infiltrator v1.8.3
  [2b0e0bc5] LanguageServer v4.5.1
  [d3d80556] LineSearches v7.2.0
⌃ [33e6dc65] MKL v0.6.3
  [687ffad2] MacroModelling v0.1.37
  [429524aa] Optim v1.9.4
  [91a5bcdd] Plots v1.40.5
  [f3b207a7] StatsPlots v0.15.7
⌅ [cf896787] SymbolServer v7.4.0

Working code

using LinearAlgebra
using MacroModelling
import StatsPlots
using Plots
using Optim, LineSearches
using Debugger, Infiltrator

results = Dict{String, Any}()

nimpstep  = 20      # Number of periods in impulse response
nplot     = 600      # Number of periods to be plotted

# introduces countercyclical markups in capital producing sector
# and variable capacity utilization through the effect on depreciation
# ENDOGENOUS EMBODIED TECHNOLOGICAL GROWTH
# overhead costs that do not depends on mt
# CAPITAL UTILIZATION
# lifecycle
# MODEL SET FOR THE ENTRY COST SHOCKS

@model Baseline begin
    ## Resource constraint
    ynet[0] = (c_yc/(ynetmcons)) * c[0] + (sk_yc/(ynetmcons)) * sk[0] + (sc_yc/(ynetmcons)) * sc[0] + actualhc_yc/(ynetmcons) * hc[0] + (ac_zc/(1-ac_zc)*actualhc_yc/(ynetmcons)) * zc[0] - (inv(zc_ac-1)*actualhc_yc/(ynetmcons)) * ac[0] + (actualhk_yc/(ynetmcons)) * hk[0] + (ak_zk/(1-ak_zk)*actualhk_yc/(ynetmcons)) * zk[0] - (inv(zk_ak-1)*actualhk_yc/(ynetmcons)) * ak[0]
    # ## Euler equation
    c[0] = c[1] - r[0]
    # ## Aggregate production function FIXME 0 or ss?
    # yc[0] = (1/(1-bb)) * chiz[0] + α * k[-1] - α * l[0] + ly[0] + ((μᶜ[ss] - 1)/(1 - bb)) * Nᶜ[0] - ((bb+μᶜ[ss]*log(Nᶜ[ss]))/(1-bb)) * μᶜ[0] + α * u[0] + (bb/(1-bb)*(θ-1)) * ac[0]
    yc[0] = (1/(1-bb)) * chiz[0] + α * k[-1] - α * l[0] + ly[0] + ((μᶜ[0] - 1)/(1 - bb)) * Nᶜ[0] - ((bb+μᶜ[0]*log(Nᶜ[0]))/(1-bb)) * μᶜ[0] + α * u[0] + (bb/(1-bb)*(θ-1)) * ac[0]
    # ## demand of capital
    yc[0] = kl[0] - l[0] + ly[0] + μᶜ[0] + d[0] * (1/(1+δ/d_pk)) + pk[0] * ((δ/(d_pk+δ))) + u[0] * edu*(1/(d_pk/δ+1))
    # ## capacity choice
    yc[0] = (1 + edp) * u[0] + μᶜ[0] + kl[0] - l[0] + ly[0] + pk[0]
    # ## labor demand in final output
    yc[0] = μᶜ[0] + ly[0] + w[0]
    # ## production of new investment goods FIXME 0 or ss?
    # yk[0] = chi[0] * (1 / ( 1 - bb)) + α * kl[0] - l[0] * (α-1/(1-lc_l)) + u[0] * α - ly[0] * (lc_l/(1-lc_l)) + Nᵏ[0] * ((μᵏ[ss]-1)/(1-bb)) - ((bb+μᵏ[ss]*log(Nᵏ[ss]))/(1-bb)) * μᵏ[0] + pk[0] * (bb/(1-bb)) + ak[0] * (bb/(1-bb)*(θ-1))
    yk[0] = chi[0] * (1 / ( 1 - bb)) + α * kl[0] - l[0] * (α-1/(1-lc_l)) + u[0] * α - ly[0] * (lc_l/(1-lc_l)) + Nᵏ[0] * ((μᵏ[0]-1)/(1-bb)) - ((bb+μᵏ[0]*log(Nᵏ[0]))/(1-bb)) * μᵏ[0] + pk[0] * (bb/(1-bb)) + ak[0] * (bb/(1-bb)*(θ-1))
    # ## Real wage
    w[0] = ζ * l[0] + c[0] + μʷ[0]
    # ## Profits embodied
    prk[0] = yk[0] + pk[0] - μᵏ[0]
    # ## Profits disembodied
    prc[0] = yc[0] - μᶜ[0]
    # ## Value of an adopted innovation for embodied
    vk[0] = ak[0] * ((1-prk_vk)) + prk[0] * prk_vk + vk[1] * (1 - prk_vk) - ak[1] * ((1-prk_vk)) - r[0] * ((1-prk_vk))
    # ## Value of an adopted innovation for disembodied
    vc[0] = ac[0] * (1-prc_vc) + prc[0] * (prc_vc) + vc[1] * (1 - prc_vc) - ac[1] * ((1-prc_vc)) - r[0] * ((1-prc_vc))
    # ## Capital accumulation
    k[0] = - u[0] * (edu*δ/(1+gk)) + k[-1] * jcof + yk[0] * (1 - jcof)
    # ## Law of motion for embodied productivity
    zk[0] = zk[-1] + sk[-1]  * (ρ*(gzk+ok)/(1+gzk)) - sf[-1] * (ρ*(gzk+ok)/(1+gzk)) + chik[-1] * ((gzk+ok)/(1+gzk))
    # ## Law of motion for disembodied productivity
    zc[0] = zc[-1] + sc[-1] * (ρ*(gzc+oc)/(1+gzc)) -  sf[-1] * (ρ*(gzc+oc)/(1+gzc))
    # ## Free entry for embodied
    sk[0] * (1 - ρ) = zk[0] - sf[0] * (ρ) + jk[1] - zk[1] - r[0]
    # ## Free entry for disembodied
    sc[0] * (1 - ρ) = zc[0] - sf[0] * (ρ) + jc[1] - zc[1] - r[0]
    # ## Bellman for not adopted disemb innovation
    - jc[0] = hc[0] * (hc_jc+ϕᶜ*elc*λᶜ/R*rz*(1-zc_ac*vc_jc)) + r[0] * ((1+hc_jc)) - zc[0] * (ϕᶜ*rz*((1-λᶜ)+λᶜ*zc_ac*vc_jc)/R) + ac[1] * (ϕᶜ*λᶜ*rz*zc_ac*vc_jc/R) - vc[1] * (ϕᶜ*λᶜ*rz*zc_ac*vc_jc/R) + sf[0] * (ϕᶜ*elc*λᶜ*rz/R*(zc_ac*vc_jc-1)) + zc[1] * (ϕᶜ*rz*(1-λᶜ)/R) - jc[1] * (ϕᶜ*rz*(1-λᶜ)/R)
    # ## law of motion for adopted disembodied innvo
    ac[0] = ac[-1] * (ϕᶜ*(1-λᶜ)/(1+gzc)) + hc[-1] * (elc*λᶜ*((ϕᶜ/(1+gzc))*zc_ac-ϕᶜ/(1+gzc))) - sf[-1] * (elc*λᶜ*((ϕᶜ/(1+gzc))*zc_ac-ϕᶜ/(1+gzc))) + zc[-1] * ((1-ϕᶜ*(1-λᶜ)/(1+gzc)))
    # ## optimal investment in adoption of disemb innov
    zc[0] = sf[0] * (1+ℓᶜ) + r[0] - hc[0] * ℓᶜ - vc[1] * (1/(1-jc_vc*ac_zc)) + ac[1] * (1/(1-jc_vc*ac_zc)) + jc[1] * (1/(vc_jc*zc_ac-1)) - zc[1] * (1/(vc_jc*zc_ac-1))
    # ## Bellman for not adopted emb innovation
    - jk[0] = hk[0] * ((hk_jk+(1-ok)*elk*λᵏ/R*ra*(1-zk_ak*vk_jk))) + r[0] * ((1+hk_jk)) - zk[0] * (ϕᵏ*ra*((1-λᵏ)+λᵏ*zk_ak*vk_jk)/R) + ak[1] * (ϕᵏ*λᵏ*ra*zk_ak*vk_jk/R) - vk[1] * (ϕᵏ*λᵏ*ra*zk_ak*vk_jk/R) + sf[0] * (ϕᵏ*elk*λᵏ*ra/R*(zk_ak*vk_jk-1)) + zk[1] * (ϕᵏ*ra*(1-λᵏ)/R) - jk[1] * (ϕᵏ*ra*(1-λᵏ)/R)
    # ## law of motion for adopted embodied innvo
    ak[0] = ak[-1] * (ϕᵏ*(1-λᵏ)/(1+gzk)) + hk[-1] * (elk*λᵏ*((ϕᵏ/(1+gzk))*zk_ak-ϕᵏ/(1+gzk))) - sf[-1] * (elk*λᵏ*((ϕᵏ/(1+gzk))*zk_ak-ϕᵏ/(1+gzk))) + zk[-1] * ((1-ϕᵏ*(1-λᵏ)/(1+gzk)))
    # ## optimal investment in adoption of emb innov
    zk[0] = sf[0] * ((1+ℓᵏ)) + r[0] - hk[0] * (ℓᵏ) - vk[1] * (1/(1-jk_vk*ak_zk)) + ak[1] * (1/(1-jk_vk*ak_zk)) + jk[1] * (1/(vk_jk*zk_ak-1)) - zk[1] * (1/(vk_jk*zk_ak-1))
    # ## Arbitrage
    pk[0] = -r[0] + d[1] * ((R-1-gpk)/R) + pk[1] * ((1+gpk)/R)
    # ## entry into final goods sector
    μᶜ[0] = - yc[0] * mucof + sf[0] * mucof + Nᶜ[0] * mucof
    # ## m
    μᶜ[0] = Nᶜ[0]  * ημᶜ
    # ## entry into capital goods sector
    μᵏ[0] = - mukcof * yk[0] - pk[0] * mukcof + sf[0] * mukcof + Nᵏ[0] * mukcof
    # ## mk
    μᵏ[0] = Nᵏ[0] * ημᵏ
    # ## equivalence between klzero and jlag
    kl[0] = k[-1]
    # ## Definition of output net of total overhead costs
    ynet[0] = yc[0] * (1/(1-oc_yc)) - Nᶜ[0] * (occ_yc/(1-oc_yc)) - Nᵏ[0] * (ock_yc/(1-oc_yc)) - sf[0] * (oc_yc/(1-oc_yc))
    # ## definition of scaling factor
    sf[0] = kl[0] + ak[0] * (bb*(1-θ)) - ac[0] * (bb*(1-θ))
    # ## definition of ynetm
    ynetm[0] = ynet[0] * (1/(1-mc_yc*inv(ynet_yc)-mk_yc*inv(ynet_yc))) - yc[0] * (mc_yc/ynetmcons) + μᶜ[0] * (mc_yc/ynetmcons) - pk[0] * (mk_yc/ynetmcons) - yk[0] * (mk_yc/ynetmcons) + μᵏ[0] * (mk_yc/ynetmcons)
    # ## Definition of total value added
    yT[0] = ynetm[0] * (ynetmcons/(ynetmcons+pkyk_yc)) + pk[0] * (pkyk_yc/(ynetmcons+pkyk_yc)) + yk[0] * (pkyk_yc/(ynetmcons+pkyk_yc))
    # ## labor demand in capital goods production
    yk[0] = - pk[0] + μᵏ[0] + w[0] + l[0] * (1/(1-lc_l)) - ly[0] * (lc_l/(1-lc_l))
    # ## embodied productivity shock process
    chi[0] = ρᵡ * chi[-1] + σᵡ* eps_chi[x]
    # ## Labor augmenting technology shock process
    chiz[0] = ρᶻᵪ * chiz[-1] + σᶻᵪ * eps_chi_z[x]
    # ## Wage markup shock process
    μʷ[0] = μʷ[-1] * ρᵐʷ + σᵐʷ * eps_muw[x]
    # ## Wage markup shock process
    chik[0] = ρᵏᵪ * chik[-1] + σᵏᵪ * eps_chi_k[x]

end

@parameters Baseline begin

    β = 0.95 # discount factor (bet)
    δ = 0.08 # depreciation rate
    ζ = 1. # labor supply curvature (fi)
    α = 1/3 # k share (al)
    g_y = 0.2 * 0.7 # ss g/y ratio
    ν = 0.6 # FIXME Comment
    θ = 1 / ν # elasticity of substitution intermediate good sector (th)
    ρ = 0.9 # parameter embodied technology
    η = 0.0
    μʷ[ss] = 1.2 # ss wage markup (muw)
    μʷₛ = 1.2 # ss wage markup (muw)
    μᶜ[ss] = 1.1 # ss wage markup (muc)
    μᶜₛ = 1.1
    Nᶜ[ss] = 1 # (nc)
    Nᶜₛ = 1 # (nc)
    Nᵏ[ss] = 1 # (nk)
    Nᵏₛ = 1 # (nk)
    dμᶜ = -μᶜₛ # (dmuc)
    ημᶜ = (dμᶜ / μᶜₛ) * Nᶜₛ
    boc = (μᶜₛ - 1)/μᶜₛ
    μᵏ[ss] = 1.2
    μᵏₛ = 1.2
    ημᵏ = ημᶜ
    λᵏ = 0.1 # (lamk)
    λᶜ = 0.1 # (lamc)
    elk = 0.9
    elc = 0.9
    ℓᵏ = elk - 1 # elasticity of λ wrt H
    ℓᶜ = elc - 1 # elasticity of λ wrt H
    o = 0.03
    oz = 0.03
    oc = 0.03
    ok = 0.03
    ϕᶜ = 1 - oc # (phic)
    ϕᵏ = 1 - ok # (phic)
    bb = 0.5 # intermediate share in final output

    ## Nonstochastic steady state
    gpk    = -0.026
    gy     =  0.024
    gk     =  gy - gpk
    gzc    = (gy-α*gk)/bb*(1-bb)/(θ-1)
    gzk    = (gpk-gzc*bb*(θ-1))/(bb*(1-θ))
    gtfp   =  gy-α*gk+gzk*(α*bb*(θ-1))/(1-α*(1-bb))
    measbls = (0.014-gy+α*gk)/(gzk*(α*bb*(θ-1))/(1-α*(1-bb)))
    gv     = gy
    gvz    = gy
    R      = (1+gy)/β
    d_pk   = R-(1+gpk)                   # definition of R
    yc_pkkc = μᶜₛ/(α*(1-bb))*(d_pk+δ) # foc for k
    yk_kk   = μᵏₛ/(α*(1-bb))*(d_pk+δ) # new capital to capital in capital production sector
    yk_k   = (gk+δ)/(1+gk)             # new capital to capital ratio
    kk_k   = yk_k/yk_kk                  # share of capital in capital production.
    kc_k   = 1-kk_k
    kk_kc  = kk_k/kc_k
    lk_lc  = kk_kc
    lk_l   = lk_lc/(lk_lc+1)
    lc_l   = 1-lk_l
    pkyk_yc= kk_kc*μᵏₛ/μᶜₛ
    mk_yc  = bb*1/θ*pkyk_yc/μᵏₛ
    mc_yc  = bb*1/θ/μᶜₛ
    pkk_yc = inv(yc_pkkc)/kc_k
    pik_yc = pkk_yc*μᶜₛ/μᵏₛ              # value of total capital stock removing fluctuations in relative price of capital due to markup variations
    prk_yc   = pkyk_yc*(1-1/θ)*bb/μᵏₛ
    prc_yc  = (1-1/θ)*bb/μᶜₛ
    prk_vk   = 1-(1+gv)*ϕᵏ/((1+gzk)*R) # bellman for va
    prc_vc = 1-(1+gvz)*ϕᶜ/((1+gzc)*R) # bellman for vz
    yc_vk    = prk_vk*inv(prk_yc)
    yc_vc   = prc_vc*inv(prc_yc)
    zk_ak   = ((gzk+ok)/(λᵏ*ϕᵏ)+1)
    zc_ac   = ((gzc+oc)/(λᶜ*ϕᶜ)+1)
    ac_zc   = inv(zc_ac)
    ak_zk   = inv(zk_ak)
    ra     = (1+gy)/(1+gzk)
    rz     = (1+gy)/(1+gzc)

    jk_yc    = (inv(1-elk*ϕᵏ*λᵏ*ra/R-(1-λᵏ)*ϕᵏ*ra/R))*(1-elk)*ϕᵏ*λᵏ*ra*zk_ak/R*inv(yc_vk) # zk * jk /yc bellman for not adopted innov
    jc_yc   = (inv(1/ϕᶜ-elc*λᶜ*rz/R-(1-λᶜ)*rz/R))*(1-elc)*λᶜ*rz*zc_ac/R*inv(yc_vc) # zc*jc/yc bellman for not adopted innov

    hk_yc    = ϕᵏ*elk*λᵏ*ra/R*(inv(yc_vk)*zk_ak-jk_yc) # zk *hk/yc
    hc_yc    = ϕᶜ*elc*λᶜ*rz/R*(inv(yc_vc)*zc_ac-jc_yc) # zc *hc/yc

    sk_yc    = jk_yc*(gzk+o)*(1+gv)*inv((1+gzk)*R) # from free entry cond't
    sc_yc   = jc_yc*(gzc+oz)*(1+gvz)*inv((1+gzc)*R)

    hc_jc  = hc_yc/jc_yc
    hk_jk  = hk_yc/jk_yc

    vc_jc  = inv(yc_vc)/jc_yc
    vk_jk  = inv(yc_vk)/jk_yc

    jc_vc=inv(vc_jc)
    jk_vk=inv(vk_jk)

    bock   = boc*pkyk_yc*(μᵏₛ-1)*μᶜₛ/(μᵏₛ*(μᶜₛ-1))
    occ_yc=boc*pik_yc
    ock_yc=bock*pik_yc
    oc_yc=occ_yc+ock_yc

    c_yc = 1-oc_yc-g_y-mc_yc-mk_yc-sk_yc-sc_yc-((ϕᶜ/(1+gzc))^2-inv(zc_ac))*hc_yc-((ϕᵏ/(1+gzk))^2-inv(zk_ak))*hk_yc
    pi_yc=(μᶜₛ-1)/μᶜₛ-oc_yc
    u[ss]=.8
    uₛ = .8
    edu=α*(1-bb)*yc_pkkc/(μᶜₛ*δ) # from foc wrt utilization, edu = elasticity of depreciation with respect to capacity
    edup=0 # partial of edu wrt u
    edp=1/3 # (edu)-1+edup/(edu*u) # elasticity of del' (i.e. elasticity of delta prima)
    actualhk_yc=hk_yc*(1-ak_zk) # total expenses in adoption of capital specific innovations
    actualhc_yc=hc_yc*(1-ac_zc) # total expenses in adoption of consumption specific innovations
    inv_Y=pkyk_yc/(pkyk_yc+1-mc_yc-mk_yc-occ_yc-ock_yc) # investment output ratio
    Y_yc=pkyk_yc/inv_Y
    ynet_yc=1-oc_yc

    # NOTE Coefficients for the log-linearization
    qcof = (1-δ)*(1+gpk)/R
    jcof = (1-δ)/(1+gk)
    vcof = (1+gy)/((1+gzk)*R)
    vzcof= (1+gy)/((1+gzc)*R)
    mucof= μᶜₛ-1
    mukcof=μᵏₛ-1

    ycof1=ynet_yc*inv(ynet_yc-mc_yc-mk_yc+pkyk_yc-actualhk_yc-actualhc_yc-sk_yc-sc_yc)
    ycof2=mc_yc*inv(ynet_yc-mc_yc-mk_yc+pkyk_yc-actualhk_yc-actualhc_yc-sk_yc-sc_yc)
    ycof3=mk_yc*inv(ynet_yc-mc_yc-mk_yc+pkyk_yc-actualhk_yc-actualhc_yc-sk_yc-sc_yc)
    ycof4=pkyk_yc*inv(ynet_yc-mc_yc-mk_yc+pkyk_yc-actualhk_yc-actualhc_yc-sk_yc-sc_yc)
    ycof=inv(ynet_yc-mc_yc-mk_yc+pkyk_yc-actualhk_yc-actualhc_yc-sk_yc-sc_yc)

    ynetmcons=1-oc_yc-mc_yc-mk_yc # fraction of ynetm in y

    # NOTE Shock to Embodied Technology
    ρᵡ = (0.7)^4   # autoregressive component
    σᵡ = 0.01    # standard deviation

    # Disembodied Technology Shock
    ρᶻᵪ = (0.7)^4
    σᶻᵪ = 0.01

    # Wage markup shock
    ρᵐʷ   = 0.60 # (rhow)
    σᵐʷ = 0.01

    # Chik shock
    ρᵏᵪ = 0.8
    σᵏᵪ = 0.01

end

model = Baseline

plot_irf(
    model,
    periods = nimpstep,
    show_plots = false,
    shocks = :all,
    save_plots = true,
    save_plots_path = "./graphs"
)

Matlab code for reference

Modifications

This is just part of the code that I want to replicate. What I did so far is to replace some variables/parameters with their unicode character when it makes sense (non-necessary). Taking the Euler equation as an example, what I did is:

% Matlab version
% Euler equation 
% 0 = c(t) - c(t+1) + r(t)
cof(2,czero) = 1;
cof(2,clead) = -1;
cof(2,rzero) = 1;

# julia version
    c[0] = c[1] - r[0]

PS: When I applied the same logic for their RBC model, I was able to replicate all the results.

Remaining part of the code

% Revision 3: This version corrects for the typo in the law of motion for adopted innovations 
% Diffusion with a two period flat lag for A, and for Z R&D and variable externality in diffusion
% Wage Markup shocks and similar production in K and Y
% growth in embodied.
% introduces countercyclical markups in capital producing sector 
% and variable capacity utilization through the effect on depreciation 
% ENDOGENOUS EMBODIED TECHNOLOGICAL GROWTH
% overhead costs that do not depends on mt
% CAPITAL UTILIZATION
% lifecycle
% MODEL SET FOR THE ENTRY COST SHOCKS

global mode

plotflag = 0;       % Set to one if you want to see plots 
                    % Hit escape to get out of a plot 
nimpstep = 500;     % number of steps in impulse response

% Model Parameters 

bet    = 0.95;       % discount factor
del    = 0.08;       % depreciation rate
fi     = 1;          % labor supply curvature
al     = 1/3;        % k share
g_y    = 0.2*0.7;    % ss g/y ratio
th     = 1/0.6;      % elasticity of substitution intermediate good sector
rho    = 0.9;        % parameter embodied technology        
eta    = 0.0;
theta  = 1/0.6;
muw    = 1.2;        % ss wage markup  
muc    = 1.1;
nc     = 1;
nk     = 1;
dmuc   = -muc;
etamuc = dmuc*nc/muc; 
boc    = (muc-1)/muc;
muk    = 1.2;
etamuk = etamuc;
lamk   = 0.1;
lamc   = 0.1;
elk    = 0.9;
elc    = 0.9;
ellk   = elk-1;
ellc   = elc-1;
o  = 0.03;
oz = 0.03;
oc = 0.03;
ok = 0.03;
phic   = 1-oc;
phik   = 1-ok;

bb     = 0.5; % intermediate share in final output

% Nonstochastic steady state

gpk    = -0.026;
gy     =  0.024;
gk     =  gy - gpk;
gzc    = (gy-al*gk)/bb*(1-bb)/(theta-1);
gzk    = (gpk-gzc*bb*(theta-1))/(bb*(1-th));
gtfp   =  gy-al*gk+gzk*(al*bb*(th-1))/(1-al*(1-bb));
measbls = (0.014-gy+al*gk)/(gzk*(al*bb*(th-1))/(1-al*(1-bb)));
gv     = gy;
gvz    = gy;
R      = (1+gy)/bet;
d_pk   = R-(1+gpk);                   % definition of R
yc_pkkc = muc/(al*(1-bb))*(d_pk+del); % foc for k
yk_kk   = muk/(al*(1-bb))*(d_pk+del); % new capital to capital in capital production sector
yk_k   = (gk+del)/(1+gk);             % new capital to capital ratio
kk_k   = yk_k/yk_kk;                  % share of capital in capital production.
kc_k   = 1-kk_k;
kk_kc  = kk_k/kc_k;
lk_lc  = kk_kc;
lk_l   = lk_lc/(lk_lc+1);
lc_l   = 1-lk_l;
pkyk_yc= kk_kc*muk/muc;
mk_yc  = bb*1/th*pkyk_yc/muk;
mc_yc  = bb*1/theta/muc;
pkk_yc = inv(yc_pkkc)/kc_k;
pik_yc = pkk_yc*muc/muk;              % value of total capital stock removing fluctuations in relative price of capital due to markup variations
prk_yc   = pkyk_yc*(1-1/th)*bb/muk;
prc_yc  = (1-1/theta)*bb/muc;
prk_vk   = 1-(1+gv)*phik/((1+gzk)*R); % bellman for va
prc_vc = 1-(1+gvz)*phic/((1+gzc)*R); % bellman for vz
yc_vk    = prk_vk*inv(prk_yc);
yc_vc   = prc_vc*inv(prc_yc);   
zk_ak   = ((gzk+ok)/(lamk*phik)+1);
zc_ac   = ((gzc+oc)/(lamc*phic)+1);
ac_zc   = inv(zc_ac);
ak_zk   = inv(zk_ak);
ra     = (1+gy)/(1+gzk);
rz     = (1+gy)/(1+gzc); 

jk_yc    = inv(1-elk*phik*lamk*ra/R-(1-lamk)*phik*ra/R)*(1-elk)*phik*lamk*ra*zk_ak/R*inv(yc_vk); % zk * jk /yc bellman for not adopted innov
jc_yc   = inv(1/phic-elc*lamc*rz/R-(1-lamc)*rz/R)*(1-elc)*lamc*rz*zc_ac/R*inv(yc_vc); % zc*jc/yc bellman for not adopted innov

hk_yc    = phik*elk*lamk*ra/R*(inv(yc_vk)*zk_ak-jk_yc); % zk *hk/yc
hc_yc    = phic*elc*lamc*rz/R*(inv(yc_vc)*zc_ac-jc_yc); % zc *hc/yc

sk_yc    = jk_yc*(gzk+o)*(1+gv)*inv((1+gzk)*R); % from free entry cond't
sc_yc   = jc_yc*(gzc+oz)*(1+gvz)*inv((1+gzc)*R);

hc_jc  = hc_yc/jc_yc; 
hk_jk  = hk_yc/jk_yc;

vc_jc  = inv(yc_vc)/jc_yc;
vk_jk  = inv(yc_vk)/jk_yc;

jc_vc=inv(vc_jc);
jk_vk=inv(vk_jk);

bock   = boc*pkyk_yc*(muk-1)*muc/(muk*(muc-1));
occ_yc=boc*pik_yc;   
ock_yc=bock*pik_yc;
oc_yc=occ_yc+ock_yc;

c_yc = 1-oc_yc-g_y-mc_yc-mk_yc-sk_yc-sc_yc-((phic/(1+gzc))^2-inv(zc_ac))*hc_yc-((phik/(1+gzk))^2-inv(zk_ak))*hk_yc;
pi_yc=(muc-1)/muc-oc_yc;        
u=.8;
edu=al*(1-bb)*yc_pkkc/(muc*del); % from foc wrt utilization, edu = elasticity of depreciation with respect to capacity
edup=0; % partial of edu wrt u
edp=1/3;%(edu)-1+edup/(edu*u); % elasticity of del' (i.e. elasticity of delta prima)
actualhk_yc=hk_yc*(1-ak_zk); % total expenses in adoption of capital specific innovations
actualhc_yc=hc_yc*(1-ac_zc); % total expenses in adoption of consumption specific innovations
inv_Y=pkyk_yc/(pkyk_yc+1-mc_yc-mk_yc-occ_yc-ock_yc); % investment output ratio
Y_yc=pkyk_yc/inv_Y;
ynet_yc=1-oc_yc;

% Coefficients for the log-linearization
qcof = (1-del)*(1+gpk)/R;
jcof = (1-del)/(1+gk);
vcof = (1+gy)/((1+gzk)*R);
vzcof= (1+gy)/((1+gzc)*R);
mucof= muc-1;
mukcof=muk-1;

ycof1=ynet_yc*(ynet_yc-mc_yc-mk_yc+pkyk_yc-actualhk_yc-actualhc_yc-sk_yc-sc_yc)^(-1);
ycof2=mc_yc*(ynet_yc-mc_yc-mk_yc+pkyk_yc-actualhk_yc-actualhc_yc-sk_yc-sc_yc)^(-1);
ycof3=mk_yc*(ynet_yc-mc_yc-mk_yc+pkyk_yc-actualhk_yc-actualhc_yc-sk_yc-sc_yc)^(-1);
ycof4=pkyk_yc*(ynet_yc-mc_yc-mk_yc+pkyk_yc-actualhk_yc-actualhc_yc-sk_yc-sc_yc)^(-1);
ycof=(ynet_yc-mc_yc-mk_yc+pkyk_yc-actualhk_yc-actualhc_yc-sk_yc-sc_yc)^(-1);

ynetmcons=1-oc_yc-mc_yc-mk_yc; % fraction of ynetm in y

% Shock to Embodied Technology 
rhochi = (0.7)^4;   % autoregressive component
chisigma = 0.01;    % standard deviation

% Disembodied Technology Shock
rhochiz = (0.7)^4;  % autoregressive component
chizsigma = 0.01;   % standard deviation

% Wage markup shock
rhomuw   = 0.60;    % autoregressive component
muwsigma = 0.01;    % standard deviation

% Chik shock
rhochik   = 0.8;    % autoregressive component
chiksigma = 0.01;   % standard deviation

%******************************************************************
%**** Begin setup for constructing log-linear system coef matrix **
%******************************************************************

nlead = 1;  % Number of leads in system 
nlag  = 1;  % Number of lags in system 

% Coefficient indicators: 
% (note xnames is not used for anything right now except to determine number of equations "neq" )

%                 1    2   3    4    5   6    7    8   9  10   11    12   13   14  15   16   17    18    19    20   21   22    23    24    25    26    27    28    29    30   31     
xnames = strvcat('yc','c','sk','sc','r','pk','hk','d','k','l','prk','prc','w','vc','vk','zk','zc','muc','nc', 'muk','nk','u', 'kl', 'ynet','sf', 'yT', 'hc','jk', 'jc', 'ak', 'ac'...
               ,  'chi','chiz','muw', 'echi','echiz','emuw', 'ly', 'yk', 'ynetm', 'chik', 'echik');
%                   32    33    34      35      36     37     38     39     40      41       42

xnum = length(xnames); %  Number of variables in system 
neq = xnum;            %  Number of equations (same)

% Ordering contemporaneous variables: 
% These are position indicators for leads and lags in constructing the coefficient matrix 
%                                                         
% Note: All shock variables must be placed last in the system!

ycpos     =  1;
cpos      =  2;
skpos     =  3;
scpos     =  4;
rpos      =  5;
pkpos     =  6;
hkpos     =  7;
dpos      =  8;
kpos      =  9;
lpos      =  10;
prkpos    =  11;
prcpos    =  12;
wpos      =  13;
vkpos     =  14;
vcpos     =  15;
zkpos     =  16;
zcpos     =  17;
mucpos    =  18;
ncpos     =  19;
mukpos    =  20;
nkpos     =  21;
upos      =  22;
klpos     =  23;
ynetpos   =  24;
sfpos     =  25;
yTpos     =  26;
hcpos     =  27;
jkpos     =  28;
jcpos     =  29;
akpos     =  30;
acpos     =  31;
lypos     =  32;
ykpos     =  33;
ynetmpos  =  34;
chipos    =  35;
chizpos   =  36;
muwpos    =  37;
chikpos   =  38;
echipos   =  39;
echizpos  =  40;
emuwpos   =  41;
echikpos  =  42;

collag  = 0;           % Position counter for start of lag coefs  
colzero = xnum;        % Position counter for start of contemp. coefs 
                       % i.e. if there is one lag and 5 equations then the "first"
                       % contemporaneous variable is at position 6 (6 = 1*5+cpos)
collead = xnum+xnum; % Position counter for start of lead coefs 

% Indicators for contemporanous coefs for each variable: 

yczero   = colzero+ycpos;
czero    = colzero+cpos;
skzero   = colzero+skpos;
hkzero   = colzero+hkpos;
rzero    = colzero+rpos;
pkzero   = colzero+pkpos;
dzero    = colzero+dpos;
kzero    = colzero+kpos;
lzero    = colzero+lpos;
prkzero  = colzero+prkpos;
wzero    = colzero+wpos;
vkzero   = colzero+vkpos; % average value of a new capital good firm prior to the vintage adjustment
zkzero   = colzero+zkpos;
zczero   = colzero+zcpos;
sczero   = colzero+scpos;
prczero  = colzero+prcpos;
vczero   = colzero+vcpos;
muczero  = colzero+mucpos;
nczero   = colzero+ncpos;
mukzero  = colzero+mukpos;
nkzero   = colzero+nkpos;
uzero    = colzero+upos;
klzero   = colzero+klpos; % klzero(t) = j(t-1)
ynetzero = colzero+ynetpos; % Yt-total overhead costs(t)
sfzero   = colzero+sfpos;
yTzero   = colzero+yTpos; % this is the total GDP in the economy accounting for the value added created at the capital producings ector.
hczero   = colzero+hcpos;
jkzero   = colzero+jkpos;
jczero   = colzero+jcpos;
akzero   = colzero+akpos;
aczero   = colzero+acpos;
chizero  = colzero+chipos;
chizzero = colzero+chizpos;
muwzero  = colzero+muwpos;
echizero  = colzero+echipos;
echizzero = colzero+echizpos;
emuwzero  = colzero+emuwpos;
lyzero    = colzero +lypos;
ykzero    = colzero +ykpos;
ynetmzero = colzero + ynetmpos;
chikzero  = colzero + chikpos;
echikzero = colzero + echikpos;

% Indicators for lead coefficients for each variable: 

yclead    = collead+ycpos;
clead    = collead+cpos;
sklead    = collead+skpos;
hklead    = collead+hkpos;
rlead    = collead+rpos;
pklead    = collead+pkpos;
dlead   = collead+dpos;
klead    = collead+kpos;
llead    = collead+lpos;
prklead   = collead+prkpos;
wlead    = collead+wpos;
vklead    = collead+vkpos;
zklead    = collead+zkpos;
zclead    = collead+zcpos;
sclead   = collead+scpos;
prclead  = collead+prcpos;
vclead   = collead+vcpos;
muclead   = collead+mucpos;
nclead    = collead+ncpos;
muklead  = collead+mukpos;
nklead   = collead+nkpos;
ulead    = collead+upos;
kllead   = collead+klpos;
ynetlead = collead+ynetpos;
sflead= collead+sfpos;
yTlead = collead+yTpos;
hclead  = collead+hcpos;
jklead  = collead+jkpos;
jclead  = collead+jcpos;
aklead  = collead+akpos;
aclead  = collead+acpos;
chilead  = collead+chipos;
chizlead = collead+chizpos;
muwlead  = collead+muwpos;
echilead  = collead+echipos;
echizlead = collead+echizpos;
emuwlead  = collead+emuwpos;
lylead  = collead+lypos;
yklead  = collead +ykpos;
ynetmlead = collead +ynetmpos;
chiklead  = collead + chikpos;
echiklead = collead + echikpos;

% Indicators for lag coefficients for each variable: 

yclag    = collag+ycpos;
clag    = collag+cpos;
sklag    = collag+skpos;
hklag    = collag+hkpos;
rlag    = collag+rpos;
pklag    = collag+pkpos;
dlag   = collag+dpos;
klag    = collag+kpos;
llag    = collag+lpos;
prklag   = collag+prkpos;
wlag    = collag+wpos;
vklag    = collag+vkpos;
zklag    = collag+zkpos;
zclag    = collag+zcpos;
sclag   = collag+scpos;
prclag  = collag+prcpos;
vclag   = collag+vcpos;
muclag   = collag+mucpos;
nclag    = collag+ncpos;
muklag  = collag+mukpos;
nklag   = collag+nkpos;
ulag    = collag+upos;
kllag   = collag+klpos;
ynetlag = collag+ynetpos;
sflag=collag+sfpos;
yTlag=collag+yTpos;
hclag  = collag+hcpos;
jklag  = collag+jkpos;
jclag  = collag+jcpos;
aklag  = collag+akpos;
aclag  = collag+acpos;
chilag  = collag+chipos;
chizlag = collag+chizpos;
muwlag  = collag+muwpos;
echilag  = collag+echipos;
echizlag = collag+echizpos;
emuwlag  = collag+emuwpos;
lylag  = collag+lypos;
yklag  = collag +ykpos;
ynetmlag = collag +ynetmpos;
chiklag  = collag + chikpos;
echiklag = collag + echikpos;

% Now I have one vector with all of the leads, leads, etc in one column and contemporaneous

% Determine number of coefficients per equation: 
ncoef = neq*(nlag+nlead+1);

cof = zeros(neq,ncoef);             % Coef matrix --- Each row is an equation

% Fill in the coefficient matrix (cof)

% * Resource constraint
% 0 = ynet(t) -1/(1-oc_y)*c_y*c(t) - 1/(1-oc_y)*th*inv_y/muk*inv(t)...
% - g_y*g(t) - 1/(1-oc_y)*s_y*s(t) - 1/(1-oc_y)*s_z*sz(t) - 1/(1-oc_y)*(1-za_z)*xz_y*xz(t)...
% +1/(1-oc_y)*za_z*xz_y*za(t) - 1/(1-oc_y)*za_z*xz_y*z(t) - 1/(1-oc_y)*(1-aa_a)*x_y*x(t)...
%+1/(1-oc_y)*aa_a*x_y*aa(t) - 1/(1-oc_y)*aa_a*x_y*a(t)

cof(1,ynetmzero)   = 1;
cof(1,czero)   = -c_yc/(ynetmcons);
%cof(1,gzero)   = -g_y/(1-oc_y);
%cof(1,sfzero)   = -.2;
cof(1,skzero)   = -sk_yc/(ynetmcons);
cof(1,sczero)  = -sc_yc/(ynetmcons);
cof(1,hczero)  = - actualhc_yc/(ynetmcons);
cof(1,zczero)  =  -ac_zc/(1-ac_zc)*actualhc_yc/(ynetmcons);
cof(1,aczero)  = inv(zc_ac-1)*actualhc_yc/(ynetmcons);
%cof(1,zllzero)  = -inv(1-za_z)*actualxz_y/(ynetmcons);
cof(1,hkzero)  = - actualhk_yc/(ynetmcons);
cof(1,zkzero)  =  -ak_zk/(1-ak_zk)*actualhk_yc/(ynetmcons);
cof(1,akzero)  = inv(zk_ak-1)*actualhk_yc/(ynetmcons);
%cof(1,allzero)  = -inv(1-aa_a)*actualx_y/(ynetmcons);

% Euler equation 
% 0 = c(t) - c(t+1) + r(t)
cof(2,czero) = 1;
cof(2,clead) = -1;
cof(2,rzero) = 1;

% Aggregate production function
% 0 = y(t) - chiz(t) al*(j(t-1)) - (1-al)*l(t) -(1/theta-1)*za(t) - (mu-1)*m(t) - (mu*log(M))*mu(t)-al*u(t)
cof(3,yczero)   = 1;
cof(3,chizzero)   = -1/(1-bb);
cof(3,klag)    = -al;
cof(3,lzero)   =   al;
cof(3,lyzero)   = -1;
cof(3,nczero)   = -(muc-1)/(1-bb);
cof(3,muczero)  = (bb+muc*log(nc))/(1-bb);
cof(3,uzero)   = -al;
cof(3,aczero) = -bb/(1-bb)*(theta-1);
%cof(3,Hyzero)   = -bb;

% * demand of capital
% 0 = y(t) - jy(t) - mu(t) - (1/(1+del*q_pk)) * pk(t) - (1/(pk_q/del+1)) * ( q(t) + edu * u(t))

cof(4,yczero)   = 1;
cof(4,klzero)  = -1;
cof(4,lzero)  =  1;
cof(4,lyzero)  = -1;
cof(4,muczero)  = -1;
cof(4,dzero)  = - (1/(1+del/d_pk));
cof(4,pkzero)   = - (del/(d_pk+del));
cof(4,uzero)   = -edu*(1/(d_pk/del+1));

% * capacity choice
% 0 = y(t) -(1+edp)*u(t) -mu(t)  -jy(t) -q(t)

cof(5,yczero) = 1;
cof(5,uzero) = -(1+edp);
cof(5,muczero) = -1;
cof(5,klzero) = -1;
cof(5,lzero) =  1;
cof(5,lyzero) = -1;
cof(5,pkzero) = -1;

% * labor demand in final output
% 0 = y(t) -ly(t) - mu(t) -w(t)

cof(6,yczero) = 1;
cof(6,muczero) = -1;
cof(6,lyzero) = -1;
cof(6,wzero) = -1;

% * production of new investment goods
% 0 = newj(t) - chi(t) - al*(1-bbk) * (jk(t) + u(t)) -(1-al)*(1-bbk)*lk(t)- bbk * Hk(t) -(muk-1) muk(t) + (muk*log(Mk)) muk(t)

cof(7,ykzero) = 1;
cof(7,chizero)  = -1/(1-bb);
cof(7,klzero) = - al;
cof(7,lzero) = al-1/(1-lc_l);
cof(7,uzero) = - al;
cof(7,lyzero) = lc_l/(1-lc_l);
cof(7,nkzero)   = -(muk-1)/(1-bb);
cof(7,mukzero)  = (bb+muk*log(nk))/(1-bb);
cof(7,pkzero)   = -bb/(1-bb);
cof(7,akzero) = -bb/(1-bb)*(th-1);

% Real wage
% 0 = w(t) - fi*l(t) - c(t) - muw(t)
cof(8,wzero)   = 1;
cof(8,lzero)   = -fi;
cof(8,czero)   = -1;
cof(8,muwzero) = -1;

% * Profits embodied
% 0 = pr(t) - newj(t) - q(t) + muk(t)
cof(9,prkzero)   = 1;
cof(9,ykzero)  = -1;
cof(9,pkzero)   = -1;
cof(9,mukzero)   = 1;

% * Profits disembodied
% 0 = prz(t) - y(t) + mu(t)
cof(10,prczero)  = 1;
cof(10,yczero)    = -1;
cof(10,muczero)   = 1;

% Value of an adopted innovation for embodied
% 0 = v(t)  - (1-vcof)*pr(t) - vcof*(v(t+1) - a(t) + a(t-1) - r(t))
cof(11,vkzero)   = 1;
cof(11,akzero)    = -(1-prk_vk);
cof(11,prkzero)  = -prk_vk;
cof(11,vklead)   = -(1-prk_vk);
cof(11,aklead)   = (1-prk_vk);
cof(11,rzero)   = (1-prk_vk);

% Value of an adopted innovation for disembodied
% 0 = vz(t)  - (1-vzcof)*prz(t) - vzcof*(vz(t+1) - z(t) + z(t-1) - r(t))
cof(12,vczero)   = 1;
cof(12,aczero)   = -(1-prc_vc);
cof(12,prczero)  = -prc_vc;
cof(12,vclead)   = -(1-prc_vc);
cof(12,aclead)   = (1-prc_vc);
cof(12,rzero)    = (1-prc_vc);

% * Capital accumulation
% 0 = j(t) + edu*del/(1+gj)* u(t) - jcof*j(t-1) - (1-jcof)*newj(t)
cof(13,kzero)   =  1;
cof(13,uzero)   =  edu*del/(1+gk);
cof(13,klag)    = -jcof;
cof(13,ykzero)  = -(1-jcof);

% Law of motion for embodied productivity
cof(14,zkzero)  =  1;
cof(14,zklag)   = -1;
cof(14,sklag)   = -rho*(gzk+ok)/(1+gzk);
cof(14,sflag)   =  rho*(gzk+ok)/(1+gzk);
cof(14,chiklag) = -(gzk+ok)/(1+gzk);

% Law of motion for disembodied productivity
cof(15,zczero)  =  1;
cof(15,zclag)   = -1;
cof(15,sclag)   = -rho*(gzc+oc)/(1+gzc);
cof(15,sflag)   =  rho*(gzc+oc)/(1+gzc);

% Free entry for embodied
% 0 = (1-rho)*s(t) - a(t-1) + rho*(q(t)+j(t-1)) - chi(t) - wa(t+1) + r(t) 
cof(16,skzero)   =  1-rho;
cof(16,zkzero)   = -1;
cof(16,sfzero)   =  rho;
cof(16,jklead)   = -1;
cof(16,zklead)   = 1;
cof(16,rzero)    = 1;

% Free entry for disembodied
% 0 = (1-rho)*sz(t) - z(t-1) + rho*(q(t)+j(t-1)) - chiz(t) - wz(t+1) + r(t) 
cof(17,sczero)    = 1-rho;
cof(17,zczero)      = -1;
cof(17,sfzero)     = rho;
cof(17,jclead)    = -1;
cof(17,zclead)     = 1;
cof(17,rzero)     = 1;

% Bellman for not adopted disemb innovation

cof(18,jczero)  = -1;
cof(18,hczero)  = -(hc_jc+phic*elc*lamc/R*rz*(1-zc_ac*vc_jc));
cof(18,rzero)   = -(1+hc_jc);
cof(18,zczero)  =  phic*rz*((1-lamc)+lamc*zc_ac*vc_jc)/R;
cof(18,aclead)  = -phic*lamc*rz*zc_ac*vc_jc/R;
cof(18,vclead)  =  phic*lamc*rz*zc_ac*vc_jc/R;
cof(18,sfzero)  = -phic*elc*lamc*rz/R*(zc_ac*vc_jc-1);
cof(18,zclead)  = -phic*rz*(1-lamc)/R;
cof(18,jclead)  =  phic*rz*(1-lamc)/R;

% law of motion for adopted disembodied innvo

cof(19,aczero) = 1; 
cof(19,aclag)  = -phic*(1-lamc)/(1+gzc);
cof(19,hclag) = -elc*lamc*((phic/(1+gzc))*zc_ac-phic/(1+gzc)); 
cof(19,sflag) = elc*lamc*((phic/(1+gzc))*zc_ac-phic/(1+gzc)); 
cof(19,zclag)= -(1-phic*(1-lamc)/(1+gzc));

% optimal investment in adoption of disemb innov

cof(20,zczero) = 1;
cof(20,sfzero) = -(1+ellc);
cof(20,rzero) = -1;
cof(20,hczero) = ellc;
cof(20,vclead) = 1/(1-jc_vc*ac_zc);
cof(20,aclead) = -1/(1-jc_vc*ac_zc);
cof(20,jclead) = -1/(vc_jc*zc_ac-1);
cof(20,zclead) =1/(vc_jc*zc_ac-1);

% Bellman for not adopted emb innovation

cof(21,jkzero)  = -1;
cof(21,hkzero)  = -(hk_jk+(1-ok)*elk*lamk/R*ra*(1-zk_ak*vk_jk));
cof(21,rzero)  = -(1+hk_jk);
cof(21,zkzero)  = phik*ra*((1-lamk)+lamk*zk_ak*vk_jk)/R;
cof(21,aklead)  = -phik*lamk*ra*zk_ak*vk_jk/R;
cof(21,vklead)  = phik*lamk*ra*zk_ak*vk_jk/R;
cof(21,sfzero)  = - phik*elk*lamk*ra/R*(zk_ak*vk_jk-1);
cof(21,zklead)  = -phik*ra*(1-lamk)/R;
cof(21,jklead)  = phik*ra*(1-lamk)/R;

% law of motion for adopted embodied innvo

cof(22,akzero) = 1; 
cof(22,aklag)  = -phik*(1-lamk)/(1+gzk);
cof(22,hklag) = -elk*lamk*((phik/(1+gzk))*zk_ak-phik/(1+gzk)); 
cof(22,sflag) = elk*lamk*((phik/(1+gzk))*zk_ak-phik/(1+gzk)); 
cof(22,zklag)= -(1-phik*(1-lamk)/(1+gzk));

% optimal investment in adoption of emb innov

cof(23,zkzero) = 1;
cof(23,sfzero) = -(1+ellk);
cof(23,rzero) = -1;
cof(23,hkzero) = ellk;
cof(23,vklead) = 1/(1-jk_vk*ak_zk);
cof(23,aklead) = -1/(1-jk_vk*ak_zk);
cof(23,jklead) = -1/(vk_jk*zk_ak-1);
cof(23,zklead) =1/(vk_jk*zk_ak-1);

% Arbitrage
% 0 = q(t) + r(t) - (R-1-gq)/R*pk(t) - (1+gq)/R*q(t+1) 
cof(24,pkzero)     = 1;
cof(24,rzero)     = 1;
cof(24,dlead)    = - (R-1-gpk)/R;
cof(24,pklead)     = -(1+gpk)/R;

% entry into final goods sector 
% 0 = mu(t) + mucof*(yf(t)+(kapa-1)*y(t)-(a+1-kapa)*m(t)-kapa*(q(t)+kl(t));+eps(t)
cof(25,muczero)    = 1;
cof(25,yczero)    = mucof;
cof(25,sfzero)    = -mucof;
cof(25,nczero)      =-mucof;

% m
% 0 = mu(t) -etamu*m(t);
cof(26,muczero)=1;
cof(26,nczero)=-etamuc;

% entry into capital goods sector 
% 0 = muk(t) + mukcof*(inv-(a+1)*mk(t)-kapa*(q(t)+kl(t));+eps(t)

cof(27,mukzero)    = 1;
cof(27,ykzero)    = mukcof;
cof(27,pkzero)    = mukcof;
cof(27,sfzero)    = -mukcof;
cof(27,nkzero)      =-mukcof;
 %cof(29,epszero)  =-mucof;

% mk
% 0 = muk(t) - etamuk * mk(t);

cof(28, mukzero) = 1;
cof(28, nkzero)  = -etamuk;

% equivalence between klzero and jlag

cof(29,klzero)=1;
cof(29,klag)=-1;

% Definition of output net of total overhead costs
% 0 = ynet(t) -1/(1-oc_y)*y(t)  + oc_y/(1-oc_y)*capsi(t)+ a/(1-exp(-a))*ocf_y)/(1-oc_y)*m(t) + a/(1-exp(-a))*ock_y)/(1-oc_y)*mk(t);

cof(30,ynetzero)   = 1;
cof(30,yczero)      = -1/(1-oc_yc);
cof(30,nczero)      = occ_yc/(1-oc_yc);
cof(30,nkzero)     = ock_yc/(1-oc_yc);
cof(30,sfzero)     = oc_yc/(1-oc_yc);

% definition of scaling factor
% 0= sf(t) -q(t) - k(t-1) +(1-muk*log(mk)) muk(t)- (muk-1) mk(t)

cof(31,sfzero) = 1;
%cof(43, qzero) = -1;
cof(31, klzero) = -1;
cof(31, akzero) = -bb*(1-th);
cof(31, aczero) = bb*(1-theta);

% definition of ynetm
% 0 = ynetm(t) 

cof(32,ynetmzero) = 1;
cof(32,ynetzero) = - 1/(1-mc_yc*inv(ynet_yc)-mk_yc*inv(ynet_yc));
cof(32,yczero) = mc_yc/ynetmcons;
cof(32,muczero) = -mc_yc/ynetmcons;
cof(32,pkzero) = mk_yc/ynetmcons;
cof(32,ykzero) = mk_yc/ynetmcons;
cof(32,mukzero) = -mk_yc/ynetmcons;

% * Definition of total value added
% 0 = yT(t) - ycof1 * ynet(t) - ycof2 * inv(t) - ycof3* muk(t)

cof(33,yTzero) =1;
cof(33,ynetmzero) =-ynetmcons/(ynetmcons+pkyk_yc);
cof(33,pkzero)  = -pkyk_yc/(ynetmcons+pkyk_yc);
cof(33,ykzero)= -pkyk_yc/(ynetmcons+pkyk_yc);

% labor demand in capital goods production
% 0 = newj(t) +q(t) -muk(t) - w(t) - 1/(1-ly_l) * l(t) + ly_l/(1-ly_l) * ly(t);
cof(34,ykzero) = 1;
cof(34,pkzero) = 1;
cof(34,mukzero) = -1;
cof(34,wzero) = -1;
cof(34,lzero) = - 1/(1-lc_l);
cof(34,lyzero) = lc_l/(1-lc_l);

% embodied productivity shock process
% 0 = chi(t)-rhochi*chi(t-1)-echi(t)
cof(35,chizero)   = 1;
cof(35,chilag)    = -rhochi;
cof(35,echizero)  = -1;

% Labor augmenting technology shock process
% 0 = chiz(t)-rhochiz*chiz(t-1)-echiz(t)
cof(36,chizzero)   =  1;
cof(36,chizlag)    = -rhochiz;
cof(36,echizzero)  = -1;

% Wage markup shock process
% 0 = muw(t)-rhomuw*muw(t-1)-emuw(t)
cof(37,muwzero)   =  1;
cof(37,muwlag)    = -rhomuw;
cof(37,emuwzero)  = -1;

% Wage markup shock process
% 0 = muw(t)-rhomuw*muw(t-1)-emuw(t)
cof(38,chikzero)   =  1;
cof(38,chiklag)    = -rhochik;
cof(38,echikzero)  = -1;

% 0 = echi(t)
cof(39,echizero)   = 1;

% 0 = echiz(t)
cof(40,echizzero)  = 1;

% 0 = emuw(t)
cof(41,emuwzero)   = 1;

% 0 = emuw(t)
cof(42,echikzero)  = 1;

%*************Begin Solution Algorithm ***************************
%*****************************************************************
%                                                                 
%  Solve a linear perfect foresight model using the gauss eig     
%  function to find the invariant subspace associated with the big 
%  roots.  This procedure will fail if the companion matrix is     
%  defective and does not have a linearly independent set of       
%  eigenvectors associated with the big roots.                     
%                                                                  
%  Input arguments:                                                
%                                                                  
%    h         Structural coef matrix (neq,neq*(nlag+1+nlag)).     
%    neq       Number of equations.                                
%    nlag      Number of lags.                                     
%    nlag      Number of lags.     (check!)!!!!!!!!!!!!!!!!                                 
%    condn     lag tolerance used as a condition number test       
%              by numeric_shift and reduced_form.                  
%    uprbnd    Inclusive upper bound for the modulus of roots      
%              allowed in the reduced form.                        
%                                                                  
%  Output arguments:                                               
%                                                                  
%    cofb      Reduced form coefficient matrix (neq,neq*nlag).     
%    scof      Observable Structure                                
%    amat      Companion form matrix                               
%    b         Contemporaneous coefficient matrix                  
%    Model satisfies:                                              
%                                                                  
%    z(t) = amat*z(t-1) + b*e(t)                                   
%                                                                  
%    where the first neq elements of z(t) are the contemporaneous  
%    values of the variables in the model                          
%    and e(t) is the shock vector of conformable dimension         
%                                                                  
%    rts       Roots returned by eig.                              
%    ia        Dimension of companion matrix (number of non-trivial
%              elements in rts).                                   
%    nexact    Number of exact shiftrights.                        
%    nnumeric  Number of numeric shiftrights.                      
%    lgroots   Number of roots greater in modulus than uprbnd.     
%    mcode     Return code: see function aimerr.                   
%*****************************************************************

% Use AIM procedure to solve model: 
uprbnd = 1+1e-8;    % Tolerance values for AIM program 
condn = 1e-8;

% ---------------------------------------------------------------------
% Run AIM
% ---------------------------------------------------------------------

[cofb,rts,ia,nex,nnum,lgrts,mcode] = ...
       aim_eig(cof,neq,nlag,nlead,condn,uprbnd);

scof = obstruct(cof,cofb,neq,nlag,nlead);

% need to calculate amat and b
% ===============================
s0 = scof(:,(neq*nlag+1):neq*(nlag+1));     %Contemp. coefs from obs. structure
amat=zeros(neq*nlag,neq*nlag);              % Initialize A matrix 
bmat=cofb(1:neq,((nlag-1)*neq+1):nlag*neq);     % Lag 1 coefficients 
i=2;
while i<=nlag;
  bmat=[bmat cofb(1:neq,((nlag-i)*neq+1):(nlag-i+1)*neq)]; % Lag i coefs 
  i=i+1;
end;
amat(1:neq,:)=bmat;                 % Coefs for equations 
if nlag>1;
 amat((length(cofb(:,1))+1):length(amat(:,1)),1:neq*(nlag-1))=eye(neq*(nlag-1));
end;
b = zeros(length(amat(:,1)),length(s0(1,:)));
b(1:length(s0(:,1)),1:length(s0(1,:))) = inv(s0);  % Store coefs 
%b=b(:,shockvec);

% COMPUTE THE COVARIANCE MATRIX FOR THE MODEL -- THEN GET CORRELATIONS 
shockvar = zeros(xnum,xnum);
shockvar(echipos,echipos)   = chisigma^2;       % variance of the shock 
shockvar(echizpos,echizpos) = chizsigma^2;      % variance of the shock 
%shockvar(egpos,egpos)       = gsigma^2;
%shockvar(eepspos,eepspos)   = epssigma^2;
shockvar(emuwpos,emuwpos)   = muwsigma^2;
omega = yvarmod(amat,b,shockvar);       % the variance covariance matrix.

domega = diag(omega)+1e-30;         % a vector of variances
sdo = domega.^0.5;
denom = sdo*sdo';
omcorr = omega./denom;              % correlation coefficients
sdratio = sdo/sdo(ycpos);
sdo = 100*sdo;

%disp(['output standard deviation (deviations)         = ' num2str(sdo(ypos))]);
%disp(['consumption standard deviation (%)             = ' num2str(sdo(cpos))]);
%disp(['investment standard deviation (%)              = ' num2str(sdo(invpos))]);
%disp(['consumption/output std dev ratio               = ' num2str(sdratio(cpos))]);
%disp(['investment/output std dev ratio                = ' num2str(sdratio(invpos))]);

% ========================================================
% Compute Impulse response function using companion form     
% solution matrix amat and contemporaneous matrix b obtained 
% from aim_run procedure: 
% ========================================================

% Shock to Embodied Technology

shock = zeros(neq,1);                        % Shock vector @
shock(echipos,1) = 1;                          % Shock variable, size of shock @
impchi = impf(amat,b,shock,nimpstep,neq);   % Call impulse respons proc @
nplot = 20;                                  % periods to be plotted
dat=(1:nplot)';                              % Date variable for plotting 

if shock(echipos,1)~=0;
if plotflag==1;

figure(1);

tfpchi=impchi(1:nplot,yTpos)-(1-al)*impchi(1:nplot,lpos)-al*impchi(1:nplot,klpos)+(1-measbls)*al*bb*(th-1)/(1-al*(1-bb))*impchi(1:nplot,akpos);
himpchi=hc_yc/(hc_yc+hk_yc)*(impchi(1:nplot,hcpos)-ac_zc/(1-ac_zc)*(impchi(1:nplot,acpos)-impchi(1:nplot,zcpos)))+ hk_yc/(hc_yc+hk_yc)*(impchi(1:nplot,hkpos)-ak_zk/(1-ak_zk)*(impchi(1:nplot,akpos)-impchi(1:nplot,zkpos)));

subplot(3,3,1);plot(dat,impchi(1:nplot,yTpos),'-');legend('Y');
subplot(3,3,2);plot(dat,impchi(1:nplot,ykpos),'-');legend('Real I');
subplot(3,3,3);plot(dat,impchi(1:nplot,lpos),'-');legend('L');
subplot(3,3,4);plot(dat,impchi(1:nplot,pkpos),'-');legend('P^K');
subplot(3,3,5);plot(dat,impchi(1:nplot,yTpos)-impchi(1:nplot,lpos),'-');legend('Y/L');
subplot(3,3,6);plot(dat,tfpchi,'-');legend('TFP');
subplot(3,3,7);plot(dat,sk_yc/(sk_yc+sc_yc)*impchi(1:nplot,skpos)+sc_yc/(sk_yc+sc_yc)*impchi(1:nplot,scpos),'-');legend('S');
subplot(3,3,8);plot(dat,himpchi,'-');legend('h');
subplot(3,3,9);plot(dat,impchi(1:nplot,upos),'-');legend('U');

end;
end;

[thorek1]

Thanks for sharing.

I can confirm the long runtime but do not have the time right now to investigate.

In the meantime you can try the following:

• only assign values to parameters (no equations referring to other parameters)
• remove the equations referring to steady state values of endogenous variables
• remove all parameter in the parameters block that are not part of the model

That should make it work fast.

Once I have a bit more time on my hand I will figure out why the code you posted worked and why it is so slow.

I hope this helps.

[gpetrini]

I tried to follow your suggestions, but still no success with both the replication and performance of the script. I also abandoned the rewrite of the code in terms of lhs = rhs and used the authors' notation (lhs - rhs = 0). I am also preserving the same parameter names as the authors to avoid that the previous results were caused by typos. Same results as before.

The rewritten code is:

using LinearAlgebra
using MacroModelling
import StatsPlots
using Plots
using Optim, LineSearches

results = Dict{String, Any}()

nimpstep  = 20      # Number of periods in impulse response
nplot     = 600      # Number of periods to be plotted

# introduces countercyclical markups in capital producing sector
# and variable capacity utilization through the effect on depreciation
# ENDOGENOUS EMBODIED TECHNOLOGICAL GROWTH
# overhead costs that do not depends on mt
# CAPITAL UTILIZATION
# lifecycle
# MODEL SET FOR THE ENTRY COST SHOCKS

@model Tmp begin
    ## Resource constraint
    ynetm[0] * 1 + c[0] * (-c_yc/(ynetmcons)) + sk[0] * ( -sk_yc/(ynetmcons) ) + sc[0] * (-sc_yc/(ynetmcons) ) + hc[0] * (- actualhc_yc/(ynetmcons) ) + zc[0] * ( -ac_zc/(1-ac_zc)*actualhc_yc/(ynetmcons) ) + ac[0] * (inv(zc_ac-1)*actualhc_yc/(ynetmcons) ) + hk[0] * (- actualhk_yc/(ynetmcons) ) + zk[0] * ( -ak_zk/(1-ak_zk)*actualhk_yc/(ynetmcons) ) + ak[0] * (inv(zk_ak-1)*actualhk_yc/(ynetmcons) ) = 0
    ## Euler equation
    c[0] * (1) + c[1] * (-1) + r[0] * (1) = 0
    ## Aggregate production function FIXME [0], [-1] or SS (muc and nc)
    yc[0] * (1) + chiz[0] * (-1/(1-bb)) + k[-1] * (-al) + l[0] * (  al) + ly[0] * (-1) + nc[0] * (-(muc[0]-1)/(1-bb)) + muc[0] * ((bb+muc[0]*log(nc[0]))/(1-bb)) + u[0] * (-al) + ac[0] * (-bb/(1-bb)*(theta-1)) = 0
    ## demand of capital
    yc[0] * (1) + kl[0] * (-1) + l[0] * ( 1) + ly[0] * (-1) + muc[0] * (-1) + d[0] * (- (1/(1+del/d_pk))) + pk[0] * (- (del/(d_pk+del))) + u[0] * (-edu*(1/(d_pk/del+1))) = 0
    ## capacity choice
    yc[0] * (1) + u[0] * (-(1+edp)) + muc[0] * (-1) + kl[0] * (-1) + l[0] * ( 1) + ly[0] * (-1) + pk[0] * (-1) = 0
    ## labor demand in final output
    yc[0] * (1) + muc[0] * (-1) + ly[0] * (-1) + w[0] * (-1) = 0
    ## production of new investment goods FIXME [0], [-1] or SS (muk and nk)
    yk[0] * (1) + chi[0]  * (-1/(1-bb)) + kl[0] * (- al) + l[0] * (al-1/(1-lc_l)) + u[0] * (- al) + ly[0] * (lc_l/(1-lc_l)) + nk[0] * (-(muk[0]-1)/(1-bb)) + muk[0] * ((bb+muk[0]*log(nk[0]))/(1-bb)) + pk[0] * (-bb/(1-bb)) + ak[0] * (-bb/(1-bb)*(th-1)) = 0
    ## Real wage
    w[0] * (1) + l[0] * (-fi) + c[0] * (-1) + muw[0] * (-1) = 0
    ## Profits embodied
    prk[0] * (1) + yk[0] * (-1) + pk[0] * (-1) + muk[0] * (1) = 0
    ## Profits disembodied
    prc[0] * (1) + yc[0] * (-1) + muc[0] * (1) = 0
    ## Value of an adopted innovation for embodied
    vk[0] * (1) + ak[0] * (-(1-prk_vk)) + prk[0] * (-prk_vk) + vk[1] * (-(1-prk_vk)) + ak[1] * ((1-prk_vk)) + r[0] * ((1-prk_vk)) = 0
    ## Value of an adopted innovation for disembodied
    vc[0] * (1) + ac[0] * (-(1-prc_vc)) + prc[0] * (-prc_vc) + vc[1] * (-(1-prc_vc)) + ac[1] * ((1-prc_vc)) + r[0] * ((1-prc_vc)) = 0
    ## Capital accumulation
    k[0] * ( 1) + u[0] * ( edu*del/(1+gk)) + k[-1] * (-jcof) + yk[0] * (-(1-jcof)) = 0
    ## Law of motion for embodied productivity
    zk[0] * ( 1) + zk[-1] * (-1) + sk[-1] * (-rho*(gzk+ok)/(1+gzk)) + sf[-1] * ( rho*(gzk+ok)/(1+gzk)) + chik[-1] * (-(gzk+ok)/(1+gzk)) = 0
    ## Law of motion for disembodied productivity
    zc[0] * ( 1) + zc[-1] * (-1) + sc[-1] * (-rho*(gzc+oc)/(1+gzc)) + sf[-1] * ( rho*(gzc+oc)/(1+gzc)) = 0
    ## Free entry for embodied
    sk[0] * ( 1-rho) + zk[0] * (-1) + sf[0] * ( rho) + jk[1] * (-1) + zk[1] * (1) + r[0] * (1) = 0
    ## Free entry for disembodied
    sc[0] * (1-rho) + zc[0] * (-1) + sf[0] * (rho) + jc[1] * (-1) + zc[1] * (1) + r[0] * (1) = 0
    ## Bellman for not adopted disemb innovation
    jc[0] * (-1) + hc[0] * (-(hc_jc+phic*elc*lamc/R*rz*(1-zc_ac*vc_jc))) + r[0] * (-(1+hc_jc)) + zc[0] * ( phic*rz*((1-lamc)+lamc*zc_ac*vc_jc)/R) + ac[1] * (-phic*lamc*rz*zc_ac*vc_jc/R) + vc[1] * ( phic*lamc*rz*zc_ac*vc_jc/R) + sf[0] * (-phic*elc*lamc*rz/R*(zc_ac*vc_jc-1)) + zc[1] * (-phic*rz*(1-lamc)/R) + jc[1] * ( phic*rz*(1-lamc)/R) = 0
    ## law of motion for adopted disembodied innvo
    ac[0] * (1) + ac[-1] * (-phic*(1-lamc)/(1+gzc)) + hc[-1] * (-elc*lamc*((phic/(1+gzc))*zc_ac-phic/(1+gzc))) + sf[-1] * (elc*lamc*((phic/(1+gzc))*zc_ac-phic/(1+gzc))) + zc[-1] * (-(1-phic*(1-lamc)/(1+gzc))) = 0
    ## optimal investment in adoption of disemb innov
    zc[0] * (1) + sf[0] * (-(1+ellc)) + r[0] * (-1) + hc[0] * (ellc) + vc[1] * (1/(1-jc_vc*ac_zc)) + ac[1] * (-1/(1-jc_vc*ac_zc)) + jc[1] * (-1/(vc_jc*zc_ac-1)) + zc[1] * (1/(vc_jc*zc_ac-1))
    ## Bellman for not adopted emb innovation
    jk[0] * (-1) + hk[0] * (-(hk_jk+(1-ok)*elk*lamk/R*ra*(1-zk_ak*vk_jk))) + r[0] * (-(1+hk_jk)) + zk[0] * (phik*ra*((1-lamk)+lamk*zk_ak*vk_jk)/R) + ak[1] * (-phik*lamk*ra*zk_ak*vk_jk/R) + vk[1] * (phik*lamk*ra*zk_ak*vk_jk/R) + sf[0] * (- phik*elk*lamk*ra/R*(zk_ak*vk_jk-1)) + zk[1] * (-phik*ra*(1-lamk)/R) + jk[1] * (phik*ra*(1-lamk)/R) = 0
    ## law of motion for adopted embodied innvo
    ak[0] * (1) + ak[-1] * (-phik*(1-lamk)/(1+gzk)) + hk[-1] * (-elk*lamk*((phik/(1+gzk))*zk_ak-phik/(1+gzk))) + sf[-1] * (elk*lamk*((phik/(1+gzk))*zk_ak-phik/(1+gzk))) + zk[-1] * (-(1-phik*(1-lamk)/(1+gzk))) = 0
    ## optimal investment in adoption of emb innov
    zk[0] * (1) + sf[0] * (-(1+ellk)) + r[0] * (-1) + hk[0] * (ellk) + vk[1] * (1/(1-jk_vk*ak_zk)) + ak[1] * (-1/(1-jk_vk*ak_zk)) + jk[1] * (-1/(vk_jk*zk_ak-1)) + zk[1] * (1/(vk_jk*zk_ak-1)) = 0
    ## Arbitrage
    pk[0] * (1) + r[0] * (1) + d[1] * (- (R-1-gpk)/R) + pk[1] * (-(1+gpk)/R) = 0
    ## entry into final goods sector
    muc[0] * (1) + yc[0] * (mucof) + sf[0] * (-mucof) + nc[0] * (-mucof) = 0
    ## m
    muc[0] * (1) + nc[0] * (-etamuc) = 0
    ## entry into capital goods sector
    muk[0] * (1) + yk[0] * (mukcof) + pk[0] * (mukcof) + sf[0] * (-mukcof) + nk[0] * (-mukcof) = 0
    ## mk
    muk[0] * (1) + nk[0] * (-etamuk) = 0
    ## equivalence between klzero and jlag
    kl[0] * (1) + k[-1] * (-1) = 0
    ## Definition of output net of total overhead costs
    ynet[0] * (1) + yc[0] * (-1/(1-oc_yc)) + nc[0] * (occ_yc/(1-oc_yc)) + nk[0] * (ock_yc/(1-oc_yc)) + sf[0] * (oc_yc/(1-oc_yc)) = 0
    ## definition of scaling factor
    sf[0] * (1) + kl[0] * (-1) + ak[0] * (-bb*(1-th)) + ac[0] * (bb*(1-theta)) = 0
    ## definition of ynetm
    ynetm[0] * (1) + ynet[0] * (- 1/(1-mc_yc*inv(ynet_yc)-mk_yc*inv(ynet_yc))) + yc[0] * (mc_yc/ynetmcons) + muc[0] * (-mc_yc/ynetmcons) + pk[0] * (mk_yc/ynetmcons) + yk[0] * (mk_yc/ynetmcons) + muk[0] * (-mk_yc/ynetmcons) = 0
    ## Definition of total value added
    yT[0] * (1) + ynetm[0] * (-ynetmcons/(ynetmcons+pkyk_yc)) + pk[0] * (-pkyk_yc/(ynetmcons+pkyk_yc)) + yk[0] * (-pkyk_yc/(ynetmcons+pkyk_yc)) = 0
    ## labor demand in capital goods production
    yk[0] * (1) + pk[0] * (1) + muk[0] * (-1) + w[0] * (-1) + l[0] * (- 1/(1-lc_l)) + ly[0] * (lc_l/(1-lc_l)) = 0
    # ## embodied productivity shock process
    chi[0] = ρᵡ * chi[-1] + σᵡ* eps_chi[x] = 0
    # ## Labor augmenting technology shock process
    chiz[0] = ρᶻᵪ * chiz[-1] + σᶻᵪ * eps_chi_z[x] = 0
    # ## Wage markup shock process
    muw[0] = muw[-1] * ρᵐʷ + σᵐʷ * eps_muw[x] = 0
    # ## Wage markup shock process
    chik[0] = ρᵏᵪ * chik[-1] + σᵏᵪ * eps_chi_k[x] = 0

end

@parameters Tmp begin
    bet    = 0.95        ## discount factor
    del    = 0.08        ## depreciation rate
    fi     = 1          ## labor supply curvature
    al     = 1/3        ## k share
    g_y    = 0.2*0.7    ## ss g/y ratio
    th     = 1/0.6      ## elasticity of substitution intermediate good sector
    rho    = 0.9        ## parameter embodied technology
    eta    = 0.0
    theta  = 1/0.6
    # muw[ss]    = 1.2        ## ss wage markup
    muw_ss    = 1.2        ## ss wage markup
    # muc[ss]    = 1.1
    muc_ss = 1.1
    # nc[ss]     = 1
    nc_ss = 1
    # nk[ss]     = 1
    nk_ss     = 1
    dmuc   = -muc_ss
    etamuc = dmuc*nc_ss/muc_ss
    boc    = (muc_ss-1)/muc_ss
    # muk[ss]    = 1.2
    muk_ss    = 1.2
    etamuk = etamuc
    lamk   = 0.1
    lamc   = 0.1
    elk    = 0.9
    elc    = 0.9
    ellk   = elk-1
    ellc   = elc-1
    o  = 0.03
    oz = 0.03
    oc = 0.03
    ok = 0.03
    phic   = 1-oc
    phik   = 1-ok
    bb     = 0.5 ## intermediate share in final output
    ## Nonstochastic steady state
    gpk    = -0.026
    gy     =  0.024
    gk     =  gy - gpk
    gzc    = (gy-al*gk)/bb*(1-bb)/(theta-1)
    gzk    = (gpk-gzc*bb*(theta-1))/(bb*(1-th))
    gtfp   =  gy-al*gk+gzk*(al*bb*(th-1))/(1-al*(1-bb))
    measbls = (0.014-gy+al*gk)/(gzk*(al*bb*(th-1))/(1-al*(1-bb)))
    gv     = gy
    gvz    = gy
    R      = (1+gy)/bet
    d_pk   = R-(1+gpk)                   ## definition of R
    yc_pkkc = muc_ss/(al*(1-bb))*(d_pk+del) ## foc for k
    yk_kk   = muk_ss/(al*(1-bb))*(d_pk+del) ## new capital to capital in capital production sector
    yk_k   = (gk+del)/(1+gk)             ## new capital to capital ratio
    kk_k   = yk_k/yk_kk                  ## share of capital in capital production.
    kc_k   = 1-kk_k
    kk_kc  = kk_k/kc_k
    lk_lc  = kk_kc
    lk_l   = lk_lc/(lk_lc+1)
    lc_l   = 1-lk_l
    pkyk_yc= kk_kc*muk_ss/muc_ss
    mk_yc  = bb*1/th*pkyk_yc/muk_ss
    mc_yc  = bb*1/theta/muc_ss
    pkk_yc = inv(yc_pkkc)/kc_k
    pik_yc = pkk_yc*muc_ss/muk_ss              ## value of total capital stock removing fluctuations in relative price of capital due to markup variations
    prk_yc   = pkyk_yc*(1-1/th)*bb/muk_ss
    prc_yc  = (1-1/theta)*bb/muc_ss
    prk_vk   = 1-(1+gv)*phik/((1+gzk)*R) ## bellman for va
    prc_vc = 1-(1+gvz)*phic/((1+gzc)*R) ## bellman for vz
    yc_vk    = prk_vk*inv(prk_yc)
    yc_vc   = prc_vc*inv(prc_yc)
    zk_ak   = ((gzk+ok)/(lamk*phik)+1)
    zc_ac   = ((gzc+oc)/(lamc*phic)+1)
    ac_zc   = inv(zc_ac)
    ak_zk   = inv(zk_ak)
    ra     = (1+gy)/(1+gzk)
    rz     = (1+gy)/(1+gzc)
    jk_yc    = inv(1-elk*phik*lamk*ra/R-(1-lamk)*phik*ra/R)*(1-elk)*phik*lamk*ra*zk_ak/R*inv(yc_vk) ## zk * jk /yc bellman for not adopted innov
    jc_yc   = inv(1/phic-elc*lamc*rz/R-(1-lamc)*rz/R)*(1-elc)*lamc*rz*zc_ac/R*inv(yc_vc) ## zc*jc/yc bellman for not adopted innov
    hk_yc    = phik*elk*lamk*ra/R*(inv(yc_vk)*zk_ak-jk_yc) ## zk *hk/yc
    hc_yc    = phic*elc*lamc*rz/R*(inv(yc_vc)*zc_ac-jc_yc) ## zc *hc/yc
    sk_yc    = jk_yc*(gzk+o)*(1+gv)*inv((1+gzk)*R) ## from free entry cond't
    sc_yc   = jc_yc*(gzc+oz)*(1+gvz)*inv((1+gzc)*R)
    hc_jc  = hc_yc/jc_yc
    hk_jk  = hk_yc/jk_yc
    vc_jc  = inv(yc_vc)/jc_yc
    vk_jk  = inv(yc_vk)/jk_yc
    jc_vc=inv(vc_jc)
    jk_vk=inv(vk_jk)
    bock   = boc*pkyk_yc*(muk_ss-1)*muc_ss/(muk_ss*(muc_ss-1))
    occ_yc=boc*pik_yc
    ock_yc=bock*pik_yc
    oc_yc=occ_yc+ock_yc
    c_yc = 1-oc_yc-g_y-mc_yc-mk_yc-sk_yc-sc_yc-((phic/(1+gzc))^2-inv(zc_ac))*hc_yc-((phik/(1+gzk))^2-inv(zk_ak))*hk_yc
    pi_yc=(muc_ss-1)/muc_ss-oc_yc
    # u[ss]=.8
    u_ss=.8
    edu=al*(1-bb)*yc_pkkc/(muc_ss*del) ## from foc wrt utilization, edu = elasticity of depreciation with respect to capacity
    edup=0 ## partial of edu wrt u
    edp=1/3##(edu)-1+edup/(edu*u); ## elasticity of del' (i.e. elasticity of delta prima)
    actualhk_yc=hk_yc*(1-ak_zk) ## total expenses in adoption of capital specific innovations
    actualhc_yc=hc_yc*(1-ac_zc) ## total expenses in adoption of consumption specific innovations
    inv_Y=pkyk_yc/(pkyk_yc+1-mc_yc-mk_yc-occ_yc-ock_yc) ## investment output ratio
    Y_yc=pkyk_yc/inv_Y
    ynet_yc=1-oc_yc
    ## Coefficients for the log-linearization
    qcof = (1-del)*(1+gpk)/R
    jcof = (1-del)/(1+gk)
    vcof = (1+gy)/((1+gzk)*R)
    vzcof= (1+gy)/((1+gzc)*R)
    mucof= muc_ss-1
    mukcof=muk_ss-1
    ycof=(ynet_yc-mc_yc-mk_yc+pkyk_yc-actualhk_yc-actualhc_yc-sk_yc-sc_yc)^(-1)
    ynetmcons=1-oc_yc-mc_yc-mk_yc ## fraction of ynetm in y
    # NOTE Shock to Embodied Technology
    ρᵡ = (0.7)^4   # autoregressive component
    σᵡ = 0.01    # standard deviation
    # Disembodied Technology Shock
    ρᶻᵪ = (0.7)^4
    σᶻᵪ = 0.01
    # Wage markup shock
    ρᵐʷ   = 0.60 # (rhow)
    σᵐʷ = 0.01
    # Chik shock
    ρᵏᵪ = 0.8
    σᵏᵪ = 0.01
end

model = Tmp

plot_irf(
    model,
    periods = nimpstep,
    show_plots = false,
    shocks = :all,
    save_plots = true,
    save_plots_path = "./graphs"
)

IRF = get_irf(
    model,
    levels = false ## in terms of SS
)

My suspicious is that I am not understanding how both matlab and this package works. As far as I understood, the matlab code set the SS values and reuses it in the equations (see, for example, muc, muk, nc, and nk).

When I use the same approach in the Julia package, how can I set initial conditions for endogenous variables? Am I missing something? I indicated with FIXMEs the two equations that I am mostly uncertain about how to proceed.

Thanks once more!

[thorek1]

here is what i meant by hard coding the parameter values (I also replaced inv(x) with 1/x, didn't test to what extent that matters). this "solves" the first order derivatives issue (I will look for a different solution so that you can still write it the way you did without having to wait for so long for derivatives). there are still some issues with the equations becayse the IRFs explode. I would recommend to recheck the timing in the equations and see that they correspond to what you have in the matlab code. be aware that the timing convention in MacroModelling is end-of-period.

using MacroModelling

@model Tmp begin
    ## Resource constraint
    ynetm[0] * 1 + c[0] * (-c_yc/(ynetmcons)) + sk[0] * ( -sk_yc/(ynetmcons) ) + sc[0] * (-sc_yc/(ynetmcons) ) + hc[0] * (- actualhc_yc/(ynetmcons) ) + zc[0] * ( -ac_zc/(1-ac_zc)*actualhc_yc/(ynetmcons) ) + ac[0] * (1/(zc_ac-1)*actualhc_yc/(ynetmcons) ) + hk[0] * (- actualhk_yc/(ynetmcons) ) + zk[0] * ( -ak_zk/(1-ak_zk)*actualhk_yc/(ynetmcons) ) + ak[0] * (1/(zk_ak-1)*actualhk_yc/(ynetmcons) ) = 0
    ## Euler equation
    c[0] * (1) + c[1] * (-1) + r[0] * (1) = 0
    ## Aggregate production function FIXME [0], [-1] or SS (muc and nc)
    yc[0] * (1) + chiz[0] * (-1/(1-bb)) + k[-1] * (-al) + l[0] * (  al) + ly[0] * (-1) + nc[0] * (-(muc[0]-1)/(1-bb)) + muc[0] * ((bb+muc[0]*log(nc[0]))/(1-bb)) + u[0] * (-al) + ac[0] * (-bb/(1-bb)*(theta-1)) = 0
    ## demand of capital
    yc[0] * (1) + kl[0] * (-1) + l[0] * ( 1) + ly[0] * (-1) + muc[0] * (-1) + d[0] * (- (1/(1+del/d_pk))) + pk[0] * (- (del/(d_pk+del))) + u[0] * (-edu*(1/(d_pk/del+1))) = 0
    ## capacity choice
    yc[0] * (1) + u[0] * (-(1+edp)) + muc[0] * (-1) + kl[0] * (-1) + l[0] * ( 1) + ly[0] * (-1) + pk[0] * (-1) = 0
    ## labor demand in final output
    yc[0] * (1) + muc[0] * (-1) + ly[0] * (-1) + w[0] * (-1) = 0
    ## production of new investment goods FIXME [0], [-1] or SS (muk and nk)
    yk[0] * (1) + chi[0]  * (-1/(1-bb)) + kl[0] * (- al) + l[0] * (al-1/(1-lc_l)) + u[0] * (- al) + ly[0] * (lc_l/(1-lc_l)) + nk[0] * (-(muk[0]-1)/(1-bb)) + muk[0] * ((bb+muk[0]*log(nk[0]))/(1-bb)) + pk[0] * (-bb/(1-bb)) + ak[0] * (-bb/(1-bb)*(th-1)) = 0
    ## Real wage
    w[0] * (1) + l[0] * (-fi) + c[0] * (-1) + muw[0] * (-1) = 0
    ## Profits embodied
    prk[0] * (1) + yk[0] * (-1) + pk[0] * (-1) + muk[0] * (1) = 0
    ## Profits disembodied
    prc[0] * (1) + yc[0] * (-1) + muc[0] * (1) = 0
    ## Value of an adopted innovation for embodied
    vk[0] * (1) + ak[0] * (-(1-prk_vk)) + prk[0] * (-prk_vk) + vk[1] * (-(1-prk_vk)) + ak[1] * ((1-prk_vk)) + r[0] * ((1-prk_vk)) = 0
    ## Value of an adopted innovation for disembodied
    vc[0] * (1) + ac[0] * (-(1-prc_vc)) + prc[0] * (-prc_vc) + vc[1] * (-(1-prc_vc)) + ac[1] * ((1-prc_vc)) + r[0] * ((1-prc_vc)) = 0
    ## Capital accumulation
    k[0] * ( 1) + u[0] * ( edu*del/(1+gk)) + k[-1] * (-jcof) + yk[0] * (-(1-jcof)) = 0
    ## Law of motion for embodied productivity
    zk[0] * ( 1) + zk[-1] * (-1) + sk[-1] * (-rho*(gzk+ok)/(1+gzk)) + sf[-1] * ( rho*(gzk+ok)/(1+gzk)) + chik[-1] * (-(gzk+ok)/(1+gzk)) = 0
    ## Law of motion for disembodied productivity
    zc[0] * ( 1) + zc[-1] * (-1) + sc[-1] * (-rho*(gzc+oc)/(1+gzc)) + sf[-1] * ( rho*(gzc+oc)/(1+gzc)) = 0
    ## Free entry for embodied
    sk[0] * ( 1-rho) + zk[0] * (-1) + sf[0] * ( rho) + jk[1] * (-1) + zk[1] * (1) + r[0] * (1) = 0
    ## Free entry for disembodied
    sc[0] * (1-rho) + zc[0] * (-1) + sf[0] * (rho) + jc[1] * (-1) + zc[1] * (1) + r[0] * (1) = 0
    ## Bellman for not adopted disemb innovation
    jc[0] * (-1) + hc[0] * (-(hc_jc+phic*elc*lamc/R*rz*(1-zc_ac*vc_jc))) + r[0] * (-(1+hc_jc)) + zc[0] * ( phic*rz*((1-lamc)+lamc*zc_ac*vc_jc)/R) + ac[1] * (-phic*lamc*rz*zc_ac*vc_jc/R) + vc[1] * ( phic*lamc*rz*zc_ac*vc_jc/R) + sf[0] * (-phic*elc*lamc*rz/R*(zc_ac*vc_jc-1)) + zc[1] * (-phic*rz*(1-lamc)/R) + jc[1] * ( phic*rz*(1-lamc)/R) = 0
    ## law of motion for adopted disembodied innvo
    ac[0] * (1) + ac[-1] * (-phic*(1-lamc)/(1+gzc)) + hc[-1] * (-elc*lamc*((phic/(1+gzc))*zc_ac-phic/(1+gzc))) + sf[-1] * (elc*lamc*((phic/(1+gzc))*zc_ac-phic/(1+gzc))) + zc[-1] * (-(1-phic*(1-lamc)/(1+gzc))) = 0
    ## optimal investment in adoption of disemb innov
    zc[0] * (1) + sf[0] * (-(1+ellc)) + r[0] * (-1) + hc[0] * (ellc) + vc[1] * (1/(1-jc_vc*ac_zc)) + ac[1] * (-1/(1-jc_vc*ac_zc)) + jc[1] * (-1/(vc_jc*zc_ac-1)) + zc[1] * (1/(vc_jc*zc_ac-1))
    ## Bellman for not adopted emb innovation
    jk[0] * (-1) + hk[0] * (-(hk_jk+(1-ok)*elk*lamk/R*ra*(1-zk_ak*vk_jk))) + r[0] * (-(1+hk_jk)) + zk[0] * (phik*ra*((1-lamk)+lamk*zk_ak*vk_jk)/R) + ak[1] * (-phik*lamk*ra*zk_ak*vk_jk/R) + vk[1] * (phik*lamk*ra*zk_ak*vk_jk/R) + sf[0] * (- phik*elk*lamk*ra/R*(zk_ak*vk_jk-1)) + zk[1] * (-phik*ra*(1-lamk)/R) + jk[1] * (phik*ra*(1-lamk)/R) = 0
    ## law of motion for adopted embodied innvo
    ak[0] * (1) + ak[-1] * (-phik*(1-lamk)/(1+gzk)) + hk[-1] * (-elk*lamk*((phik/(1+gzk))*zk_ak-phik/(1+gzk))) + sf[-1] * (elk*lamk*((phik/(1+gzk))*zk_ak-phik/(1+gzk))) + zk[-1] * (-(1-phik*(1-lamk)/(1+gzk))) = 0
    ## optimal investment in adoption of emb innov
    zk[0] * (1) + sf[0] * (-(1+ellk)) + r[0] * (-1) + hk[0] * (ellk) + vk[1] * (1/(1-jk_vk*ak_zk)) + ak[1] * (-1/(1-jk_vk*ak_zk)) + jk[1] * (-1/(vk_jk*zk_ak-1)) + zk[1] * (1/(vk_jk*zk_ak-1)) = 0
    ## Arbitrage
    pk[0] * (1) + r[0] * (1) + d[1] * (- (R-1-gpk)/R) + pk[1] * (-(1+gpk)/R) = 0
    ## entry into final goods sector
    muc[0] * (1) + yc[0] * (mucof) + sf[0] * (-mucof) + nc[0] * (-mucof) = 0
    ## m
    muc[0] * (1) + nc[0] * (-etamuc) = 0
    ## entry into capital goods sector
    muk[0] * (1) + yk[0] * (mukcof) + pk[0] * (mukcof) + sf[0] * (-mukcof) + nk[0] * (-mukcof) = 0
    ## mk
    muk[0] * (1) + nk[0] * (-etamuk) = 0
    ## equivalence between klzero and jlag
    kl[0] * (1) + k[-1] * (-1) = 0
    ## Definition of output net of total overhead costs
    ynet[0] * (1) + yc[0] * (-1/(1-oc_yc)) + nc[0] * (occ_yc/(1-oc_yc)) + nk[0] * (ock_yc/(1-oc_yc)) + sf[0] * (oc_yc/(1-oc_yc)) = 0
    ## definition of scaling factor
    sf[0] * (1) + kl[0] * (-1) + ak[0] * (-bb*(1-th)) + ac[0] * (bb*(1-theta)) = 0
    ## definition of ynetm
    ynetm[0] * (1) + ynet[0] * (- 1/(1-mc_yc*1/(ynet_yc)-mk_yc*1/(ynet_yc))) + yc[0] * (mc_yc/ynetmcons) + muc[0] * (-mc_yc/ynetmcons) + pk[0] * (mk_yc/ynetmcons) + yk[0] * (mk_yc/ynetmcons) + muk[0] * (-mk_yc/ynetmcons) = 0
    ## Definition of total value added
    yT[0] * (1) + ynetm[0] * (-ynetmcons/(ynetmcons+pkyk_yc)) + pk[0] * (-pkyk_yc/(ynetmcons+pkyk_yc)) + yk[0] * (-pkyk_yc/(ynetmcons+pkyk_yc)) = 0
    ## labor demand in capital goods production
    yk[0] * (1) + pk[0] * (1) + muk[0] * (-1) + w[0] * (-1) + l[0] * (- 1/(1-lc_l)) + ly[0] * (lc_l/(1-lc_l)) = 0
    # ## embodied productivity shock process
    chi[0] = ρᵡ * chi[-1] + σᵡ* eps_chi[x] = 0
    # ## Labor augmenting technology shock process
    chiz[0] = ρᶻᵪ * chiz[-1] + σᶻᵪ * eps_chi_z[x] = 0
    # ## Wage markup shock process
    muw[0] = muw[-1] * ρᵐʷ + σᵐʷ * eps_muw[x] = 0
    # ## Wage markup shock process
    chik[0] = ρᵏᵪ * chik[-1] + σᵏᵪ * eps_chi_k[x] = 0
end

@parameters Tmp begin
    bet    = 0.95        ## discount factor
    del    = 0.08        ## depreciation rate
    fi     = 1          ## labor supply curvature
    al     = 1/3        ## k share
    g_y    = 0.2*0.7    ## ss g/y ratio
    th     = 1/0.6      ## elasticity of substitution intermediate good sector
    rho    = 0.9        ## parameter embodied technology
    eta    = 0.0
    theta  = 1/0.6
    # muw[ss]    = 1.2        ## ss wage markup
    muw_ss    = 1.2        ## ss wage markup
    # muc[ss]    = 1.1
    muc_ss = 1.1
    # nc[ss]     = 1
    nc_ss = 1
    # nk[ss]     = 1
    nk_ss     = 1
    dmuc   = -1.1 # -muc_ss
    etamuc = -1 # dmuc*nc_ss/muc_ss
    boc    = 0.09090909090909098 # (muc_ss-1)/muc_ss
    # muk[ss]    = 1.2
    muk_ss    = 1.2
    etamuk = -1 # etamuc
    lamk   = 0.1
    lamc   = 0.1
    elk    = 0.9
    elc    = 0.9
    ellk   = -0.09999999999999998 # elk-1
    ellc   = -0.09999999999999998 # elc-1
    o  = 0.03
    oz = 0.03
    oc = 0.03
    ok = 0.03
    phic   = 0.97 # 1-oc
    phik   = 0.97 # 1-ok
    bb     = 0.5 ## intermediate share in final output
    ## Nonstochastic steady state
    gpk    = -0.026
    gy     =  0.024
    gk     =  .05 # gy - gpk
    gzc    = .011 # (gy-al*gk)/bb*(1-bb)/(theta-1)
    gzk    = .089 # (gpk-gzc*bb*(theta-1))/(bb*(1-th))
    gtfp   =  .0192 # gy-al*gk+gzk*(al*bb*(th-1))/(1-al*(1-bb))
    measbls = 0.5617977528089887 # (0.014-gy+al*gk)/(gzk*(al*bb*(th-1))/(1-al*(1-bb)))
    gv     = 0.024 # gy
    gvz    = 0.024 # gy
    R      = 1.0778947368421052 # (1+gy)/bet
    d_pk   = 0.10389473684210526 # R-(1+gpk)                   ## definition of R
    yc_pkkc = 1.2137052631578948 # muc_ss/(al*(1-bb))*(d_pk+del) ## foc for k
    yk_kk   = 1.3240421052631581 # muk_ss/(al*(1-bb))*(d_pk+del) ## new capital to capital in capital production sector
    yk_k   = 0.12380952380952381 # (gk+del)/(1+gk)             ## new capital to capital ratio
    kk_k   = 0.09350875120766361 # yk_k/yk_kk                  ## share of capital in capital production.
    kc_k   = 0.9064912487923364 # 1-kk_k
    kk_kc  = 0.1031546099669905 # kk_k/kc_k
    lk_lc  = 0.1031546099669905 # kk_kc
    lk_l   = 0.09350875120766361 # lk_lc/(lk_lc+1)
    lc_l   = 0.9064912487923364 # 1-lk_l
    pkyk_yc= 0.11253230178217144 # kk_kc*muk_ss/muc_ss
    mk_yc  = 0.02813307544554286 # bb*1/th*pkyk_yc/muk_ss
    mc_yc  = 0.2727272727272727 # bb*1/theta/muc_ss
    pkk_yc = 0.9089147451636925 # 1/(yc_pkkc)/kc_k
    pik_yc = 0.833171849733385 # pkk_yc*muc_ss/muk_ss              ## value of total capital stock removing fluctuations in relative price of capital due to markup variations
    prk_yc   = 0.01875538363036191 # pkyk_yc*(1-1/th)*bb/muk_ss
    prc_yc  = 0.18181818181818182 # (1-1/theta)*bb/muc_ss
    prk_vk   = 0.15381083562901743 # 1-(1+gv)*phik/((1+gzk)*R) ## bellman for va
    prc_vc = 0.08852621167161212 # 1-(1+gvz)*phic/((1+gzc)*R) ## bellman for vz
    yc_vk    = 8.200889870363556 # prk_vk*1/(prk_yc)
    yc_vc   = 0.48689416419386666 # prc_vc*1/(prc_yc)
    zk_ak   = 2.2268041237113403 # ((gzk+ok)/(lamk*phik)+1)
    zc_ac   = 1.4226804123711339 # ((gzc+oc)/(lamc*phic)+1)
    ac_zc   = 0.7028985507246378 # 1/(zc_ac)
    ak_zk   = 0.44907407407407407 # 1/(zk_ak)
    ra     = 0.9403122130394859 # (1+gy)/(1+gzk)
    rz     = 1.0128585558852623 # (1+gy)/(1+gzc)
    jk_yc    = 0.014159338410620156 # 1/(1-elk*phik*lamk*ra/R-(1-lamk)*phik*ra/R)*(1-elk)*phik*lamk*ra*zk_ak/R*1/(yc_vk) ## zk * jk /yc bellman for not adopted innov
    jc_yc   = 0.2727626948903887 # 1/(1/phic-elc*lamc*rz/R-(1-lamc)*rz/R)*(1-elc)*lamc*rz*zc_ac/R*1/(yc_vc) ## zc*jc/yc bellman for not adopted innov
    hk_yc    = 0.019600737056023755 # phik*elk*lamk*ra/R*(1/(yc_vk)*zk_ak-jk_yc) ## zk *hk/yc
    hc_yc    = 0.21731983257587298 # phic*elc*lamc*rz/R*(1/(yc_vc)*zc_ac-jc_yc) ## zc *hc/yc
    sk_yc    = 0.0014698927523605224 # jk_yc*(gzk+o)*(1+gv)*1/((1+gzk)*R) ## from free entry cond't
    sc_yc   = 0.01050851331946651 # jc_yc*(gzc+oz)*(1+gvz)*1/((1+gzc)*R)
    hc_jc  = 0.7967359050445085 # hc_yc/jc_yc
    hk_jk  = 1.3842975206611559 # hk_yc/jk_yc
    vc_jc  = 7.529748283752845 # 1/(yc_vc)/jc_yc
    vk_jk  = 8.61184210526315 # 1/(yc_vk)/jk_yc
    jc_vc= 0.13280656435192248 # 1/(vc_jc)
    jk_vk= 0.11611917494270446 # 1/(vk_jk)
    bock   = 0.018755383630361902 # boc*pkyk_yc*(muk_ss-1)*muc_ss/(muk_ss*(muc_ss-1))
    occ_yc= 0.07574289543030778 # boc*pik_yc
    ock_yc= 0.015626457671767874 # bock*pik_yc
    oc_yc= 0.09136935310207565 # occ_yc+ock_yc
    c_yc = 0.40174590233878776 # 1-oc_yc-g_y-mc_yc-mk_yc-sk_yc-sc_yc-((phic/(1+gzc))^2-1/(zc_ac))*hc_yc-((phik/(1+gzk))^2-1/(zk_ak))*hk_yc
    pi_yc= -0.00046026219298467286 # (muc_ss-1)/muc_ss-oc_yc
    # u[ss]=.8
    u_ss=.8
    edu= 2.2986842105263157 # al*(1-bb)*yc_pkkc/(muc_ss*del) ## from foc wrt utilization, edu = elasticity of depreciation with respect to capacity
    edup=0 ## partial of edu wrt u
    edp=1/3##(edu)-1+edup/(edu*u); ## elasticity of del' (i.e. elasticity of delta prima)
    actualhk_yc= 0.010798554211420494 # hk_yc*(1-ak_zk) ## total expenses in adoption of capital specific innovations
    actualhc_yc= 0.06456603721457094 # hc_yc*(1-ac_zc) ## total expenses in adoption of consumption specific innovations
    inv_Y= 0.1562292038136742 # pkyk_yc/(pkyk_yc+1-mc_yc-mk_yc-occ_yc-ock_yc) ## investment output ratio
    Y_yc= 0.7203026005072803 # pkyk_yc/inv_Y
    ynet_yc= 0.9086306468979244 # 1-oc_yc
    ## Coefficients for the log-linearization
    qcof = 0.83132421875 # (1-del)*(1+gpk)/R
    jcof = 0.8761904761904762 # (1-del)/(1+gk)
    vcof =  0.8723599632690543 # (1+gy)/((1+gzk)*R)
    vzcof= 0.9396636993076164 # (1+gy)/((1+gzc)*R)
    mucof= 0.10000000000000009 # muc_ss-1
    mukcof= 0.19999999999999996 # muk_ss-1
    ycof= 1.5798796562140973 # (ynet_yc-mc_yc-mk_yc+pkyk_yc-actualhk_yc-actualhc_yc-sk_yc-sc_yc)^(-1)
    ynetmcons= 0.6077702987251088 # 1-oc_yc-mc_yc-mk_yc ## fraction of ynetm in y
    # NOTE Shock to Embodied Technology
    ρᵡ = (0.7)^4   # autoregressive component
    σᵡ = 0.01    # standard deviation
    # Disembodied Technology Shock
    ρᶻᵪ = (0.7)^4
    σᶻᵪ = 0.01
    # Wage markup shock
    ρᵐʷ   = 0.60 # (rhow)
    σᵐʷ = 0.01
    # Chik shock
    ρᵏᵪ = 0.8
    σᵏᵪ = 0.01
end