BASEforHANK · BenjaminBorn · Mar 14, 2024 · Mar 14, 2024 · Mar 14, 2024
diff --git a/src/Model/input_aggregate_names.jl b/src/Model/input_aggregate_names.jl
@@ -4,7 +4,7 @@
 # contained in the aggregate model code as they are parameter free but change whenever the 
 # distribution changes and do not show up in any aggregate model equation.
 
-n_rep = 0 # Number of n_rep of some model equations (e.g. countries, industries)
+n_rep = 1 # Number of n_rep of some model equations (e.g. countries, industries)
 
 # List of aggregate shocks, without duplication (e.g. across countries or industries)
 shock_names = [:Z, :ZI, :μ, :μw, :A, :Rshock, :Gshock, :Tprogshock, :Sshock]
@@ -94,20 +94,35 @@ dup_control_names = ["",
 
 state_names = setdiff(state_names,dup_state_names)
 for  j = 1:n_rep # Create names of duplicated state variables (e.g. across countries or industries)
-    append!(state_names,dup_state_names.*string(j))
+    if j==1
+        aux = ""
+    else
+        aux = string(j)
+    end
+    append!(state_names,dup_state_names.*aux)
 end
 
 shock_states = string.(shock_names)
 dup_shock_states = intersect(string.(shock_names), dup_state_names)
 shock_states = setdiff(shock_states,dup_shock_states)
 for  j = 1:n_rep # Create names of duplicated shock variables (e.g. across countries or industries)
-    append!(shock_states ,dup_shock_states .*string(j))
+    if j==1
+        aux = ""
+    else
+        aux = string(j)
+    end
+    append!(shock_states ,dup_shock_states .*aux)
 end
 shock_names = Symbol.(shock_states)
 
 control_names= setdiff(control_names,dup_control_names)
 for j = 1:n_rep # Create names of duplicated control variables (e.g. across countries or industries)
-    append!(control_names,dup_control_names.*string(j))
+    if j==1
+        aux = ""
+    else
+        aux = string(j)
+    end
+    append!(control_names,dup_control_names.*aux)
 end
 
 # All controls in one array

diff --git a/src/Model/input_aggregate_steady_state.jl b/src/Model/input_aggregate_steady_state.jl
@@ -1,72 +1,73 @@
+@R$1 #no replication 
 # Set aggregate steady state variabel values
-ASS = 1.0
-ZSS = 1.0
-ZISS = 1.0
-μSS = m_par.μ
-μwSS = m_par.μw
-τprogSS = m_par.τprog 
-τlevSS = m_par.τlev
+A$SS = 1.0
+Z$SS = 1.0
+ZI$SS = 1.0
+μ$SS = m_par.μ
+μw$SS = m_par.μw
+τprog$SS = m_par.τprog 
+τlev$SS = m_par.τlev
 
-σSS = 1.0
-τprog_obsSS = 1.0
-GshockSS = 1.0
-RshockSS = 1.0
-TprogshockSS = 1.0
+σ$SS = 1.0
+τprog_obs$SS = 1.0
+Gshock$SS = 1.0
+Rshock$SS = 1.0
+Tprogshock$SS = 1.0
 
-SshockSS = 1.0
-RBSS = m_par.RB
-LPSS = 1 + rSS - RBSS
-LPXASS = 1 + rSS - RBSS
-ISS = m_par.δ_0 * KSS
+Sshock$SS = 1.0
+RB$SS = m_par.RB
+LP$SS = 1 + rSS - RBSS
+LPX$ASS = 1 + rSS - RBSS
+I$SS = m_par.δ_0 * KSS
 
-πSS = 1.0
-πwSS = 1.0
+π$SS = 1.0
+πw$SS = 1.0
 
 BDSS = -sum(distr_mSS .* (n_par.grid_m .< 0) .* n_par.grid_m)
 # Calculate taxes and government expenditures
-TSS = dot(distr_ySS, taxrev) + av_tax_rateSS * ((1.0 .- 1.0 ./ m_par.μw) .* wSS .* NSS)
+T$SS = dot(distr_y$SS, taxrev) + av_tax_rate$SS * ((1.0 .- 1.0 ./ m_par.μw) .* w$SS .* N$SS)
 
-BgovSS = BSS .- qΠSS_fnc(YSS, m_par.RB, m_par) .+ 1.0
-GSS = TSS - (m_par.RB ./ m_par.π - 1.0) * BgovSS
-mcSS = 1.0 ./ m_par.μ
-firm_profitsSS = (1.0 - mcSS) .* YSS
-qΠSS = qΠSS_fnc(YSS, RBSS, m_par)
-qΠlagSS = qΠSS
-RLSS = m_par.RB
+Bgov$SS = B$SS .- qΠSS_fnc(Y$SS, m_par.RB, m_par) .+ 1.0
+G$SS = T$SS - (m_par.RB ./ m_par.π - 1.0) * Bgov$SS
+mc$SS = 1.0 ./ m_par.μ
+firm_profits$SS = (1.0 - mc$SS) .* Y$SS
+qΠ$SS = qΠSS_fnc(Y$SS, RB$SS, m_par)
+qΠlag$SS = qΠ$SS
+RL$SS = m_par.RB
 
-CSS = (YSS - m_par.δ_0 * KSS - GSS - m_par.Rbar * BDSS)
+C$SS = (Y$SS - m_par.δ_0 * K$SS - G$SS - m_par.Rbar * BD$SS)
 
-qSS = 1.0
-mcwSS = 1.0 ./ m_par.μw
-mcwwSS = wSS * mcwSS
-uSS = 1.0
-unionprofitsSS = (1.0 - mcwSS) .* wSS .* NSS
+q$SS = 1.0
+mcw$SS = 1.0 ./ m_par.μw
+mcww$SS = w$SS * mcw$SS
+u$SS = 1.0
+unionprofits$SS = (1.0 - mcw$SS) .* w$SS .* N$SS
 
-BYSS = BSS / YSS
-TYSS = TSS / YSS
-TlagSS = TSS
+BY$SS = B$SS / Y$SS
+TY$SS = T$SS / Y$SS
+Tlag$SS = T$SS
 
-YlagSS = YSS
-BgovlagSS = BgovSS
-GlagSS = GSS
-IlagSS = ISS
-wlagSS = wSS
-qlagSS = qSS
-ClagSS = CSS
-av_tax_ratelagSS = av_tax_rateSS
-τproglagSS = τprogSS
+Ylag$SS = Y$SS
+Bgovlag$SS = Bgov$SS
+Glag$SS = G$SS
+Ilag$SS = I$SS
+wlag$SS = w$SS
+qlag$SS = q$SS
+Clag$SS = C$SS
+av_tax_ratelag$SS = av_tax_rate$SS
+τproglag$SS = τprog$SS
 
-YgrowthSS = 1.0
-BgovgrowthSS = 1.0
-IgrowthSS = 1.0
-wgrowthSS = 1.0
-CgrowthSS = 1.0
-TgrowthSS = 1.0
-HtSS = 1.0
+Ygrowth$SS = 1.0
+Bgovgrowth$SS = 1.0
+Igrowth$SS = 1.0
+wgrowth$SS = 1.0
+Cgrowth$SS = 1.0
+Tgrowth$SS = 1.0
+Ht$SS = 1.0
 
 # calculate and display aggregate moments
-println("BSS/YSS: ", BSS / YSS)
-println("StockshareSS: ", (qΠSS_fnc(YSS, m_par.RB, m_par) .- 1.0) / BSS)
-println("BgovSS/YSS: ", BgovSS / YSS)
-println("GSS/YSS: ", GSS / YSS)
-println("firm_profitsSS: ", firm_profitsSS)
+println("B$SS/Y$SS: ", B$SS / Y$SS)
+println("Stockshare$SS: ", (qΠSS_fnc(Y$SS, m_par.RB, m_par) .- 1.0) / B$SS)
+println("Bgov$SS/Y$SS: ", Bgov$SS / Y$SS)
+println("G$SS/Y$SS: ", G$SS / Y$SS)
+println("firm_profits$SS: ", firm_profits$SS)
diff --git a/src/Preprocessor/PreprocessInputs.jl b/src/Preprocessor/PreprocessInputs.jl
@@ -12,13 +12,8 @@ insert_index = findall(x -> x == "    # aggregate model marker", model_template_
 
 model_input_file = open("Model/input_aggregate_model.jl")
 model_input_text = read(model_input_file, String)
-model_input_file = open("Model/input_aggregate_model.jl")
-model_input_lines = readlines(model_input_file)
 
 deblank(S::String) = filter(x -> !isspace(x), S)
-number_of_equations =
-    length(collect(eachmatch(r"F\[indexes.", deblank(model_input_text)))) -
-    length(collect(eachmatch(r"\#F\[indexes.", deblank(model_input_text))))
 
 open("Preprocessor/generated_fcns/FSYS_agg_generated.jl", "w") do h
     println(
@@ -31,24 +26,32 @@ open("Preprocessor/generated_fcns/FSYS_agg_generated.jl", "w") do h
     for i = 1:insert_index-1
         println(h, model_template_lines[i])
     end
+
     println(h, "\n")
-    for i in model_input_lines  # input the model text
-        if occursin("@R", i)    # if there is a repetition marker create replications
-            n1      = findfirst("@R", i)
-            rsym    = i[n1[end]+1]
-            n2      = findfirst(" ", i[n1[end]+1:end])
-            repl    = parse(Int,i[n1[end]+2:n1[end]+n2[end]])
-            for j = 1:repl
-                line = i[n1[end]+n2[end]+1:end]
-                line = replace(line, rsym => string(j)) 
-                println(h, line)
+    # produce replications according to magic comment
+    if occursin("@R",model_input_text)
+        n1      = findfirst("@R", model_input_text)
+        rsym    = model_input_text[n1[end]+1]
+        n2      = findfirst(r"[^\d]", model_input_text[n1[end]+2:end])
+        n_rep   = parse(Int,model_input_text[n1[end]+2:n1[end]+n2[end]+1])
+        for j = 1:n_rep
+            if j==1
+                aux = ""
+                write(h, model_input_text[1:n1[1]-1])
+                text = model_input_text[n1[end]+n2[end]+1:end]
+            else
+                aux = string(j)
+                text = "\n#------------------------------------------------------------------------------\n"*
+                       "# Economy/Sector " * string(j) * " starts here" *
+                       "\n#------------------------------------------------------------------------------\n"*  
+                       model_input_text[n1[end]+n2[end]+1:end]
             end
-            if occursin("F[indexes.", i)
-                number_of_equations += repl-1
-            end
-        else
-            println(h, i)
+
+            text = replace(text, rsym => aux) 
+            write(h, text)
         end
+    else
+        write(h, model_input_text)
     end
     println(h, "\n")
     for i = insert_index+1:length(model_template_lines)
@@ -58,15 +61,22 @@ end
 close(model_template_file)
 close(model_input_file)
 
+# Check number of equations
+Fsys_agg_file = open("Preprocessor/generated_fcns/FSYS_agg_generated.jl")
+Fsys_agg_lines= readlines(Fsys_agg_file)
+Fsys_agg_lines_unique = unique(Fsys_agg_lines)
+number_of_equations =
+    count(occursin.("F[indexes.", deblank.(Fsys_agg_lines_unique))) -
+    count(occursin.("#F[indexes.", deblank.(Fsys_agg_lines_unique)))
+
+
 # aggregate steady state
 SS_template_file    = open("Preprocessor/template_fcns/prepare_linearization.jl")
 SS_template_lines   = readlines(SS_template_file)
 insert_index        = findall(x -> x == "    # aggregate steady state marker", SS_template_lines)[1]
 
 SS_input_file       = open("Model/input_aggregate_steady_state.jl")
 SS_input_text       = read(SS_input_file, String)
-SS_input_file       = open("Model/input_aggregate_steady_state.jl")
-SS_input_lines      = readlines(SS_input_file)
 
 open("Preprocessor/generated_fcns/prepare_linearization_generated.jl", "w") do h
     println(
@@ -80,20 +90,27 @@ open("Preprocessor/generated_fcns/prepare_linearization_generated.jl", "w") do h
     end
     println(h, "\n")
     println(h, "@set! n_par.n_agg_eqn = $number_of_equations")
-    for i in SS_input_lines  # input the model text
-        if occursin("@R", i)    # if there is a repetition marker create replications
-            n1      = findfirst("@R", i)
-            rsym    = i[n1[end]+1]
-            n2      = findfirst(" ", i[n1[end]+1:end])
-            repl    = parse(Int,i[n1[end]+2:n1[end]+n2[end]])
-            for j = 1:repl
-                line = i[n1[end]+n2[end]+1:end]
-                line = replace(line, rsym => string(j)) 
-                println(h, line)
+    if occursin("@R",SS_input_text)
+        n1      = findfirst("@R", SS_input_text)
+        rsym    = SS_input_text[n1[end]+1]
+        n2      = findfirst(r"[^\d]", SS_input_text[n1[end]+2:end])
+        n_rep   = parse(Int,SS_input_text[n1[end]+2:n1[end]+n2[end]+1])
+        for j = 1:n_rep
+            if j==1
+                aux = ""
+                text = SS_input_text[1:n1[1]-1]*SS_input_text[n1[end]+n2[end]+1:end]
+            else
+                aux = string(j)
+                text = "\n#------------------------------------------------------------------------------\n"*
+                       "# Economy/Sector " * string(j) * " starts here" *
+                       "\n#------------------------------------------------------------------------------\n"*  
+                       SS_input_text[n1[end]+n2[end]+1:end]
             end
-        else
-            println(h, i)
+            text = replace(text, rsym => aux) 
+            write(h, text)
         end
+    else
+        write(h, SS_input_text)
     end
     println(h, "\n")
     for i = insert_index+1:length(SS_template_lines)

diff --git a/src/SubModules/PerturbationSolution/FSYS.jl b/src/SubModules/PerturbationSolution/FSYS.jl
@@ -47,7 +47,6 @@ function Fsys(
     ############################################################################
     #            I. Read out argument values                                   #
     ############################################################################
-    # rougly 10% of computing time, more if uncompress is actually called
 
     ############################################################################
     # I.1. Generate code that reads aggregate states/controls
@@ -59,7 +58,6 @@ function Fsys(
     # @generate_equations(aggr_names)
     @generate_equations()
 
-
     ############################################################################
     # I.2. Read out  perturbed distributions
     ############################################################################
@@ -155,12 +153,12 @@ function Fsys(
 
     # Average Human Capital =
     # average productivity (at the productivit grid, used to normalize to 0)
-    tax_prog_scale = (m_par.γ + m_par.τprog ) / ((m_par.γ + τprog))
-    H = dot(distr_y[1:end-1], n_par.grid_y[1:end-1])
-    KP = dot(n_par.grid_k, distr_k[:])
-    Htact = dot(distr_y[1:end-1], (n_par.grid_y[1:end-1] / H) .^ (tax_prog_scale))
-    BP = dot(n_par.grid_m, distr_m[:])
-    BDact = -sum(distr_m .* (n_par.grid_m .< 0) .* n_par.grid_m)
+    tax_prog_scale  = (m_par.γ + m_par.τprog ) / ((m_par.γ + τprog))
+    H               = dot(distr_y[1:end-1], n_par.grid_y[1:end-1])
+    KP              = dot(n_par.grid_k, distr_k[:])
+    Htact           = dot(distr_y[1:end-1], (n_par.grid_y[1:end-1] / H) .^ (tax_prog_scale))
+    BP              = dot(n_par.grid_m, distr_m[:])
+    BDact           = -sum(distr_m .* (n_par.grid_m .< 0) .* n_par.grid_m)
 
     ############################################################################
     #               III. 2. Heterogeneous Agent Part                           #
@@ -187,19 +185,17 @@ function Fsys(
     )
 
     # Calculate Taxes
-    tax_prog_scale = (m_par.γ + m_par.τprog ) / ((m_par.γ + τprog))
-    LC = mcw * w .* N ./ Ht
+    tax_prog_scale  = (m_par.γ + m_par.τprog ) / ((m_par.γ + τprog))
+    LC              = mcw * w .* N ./ Ht
     taxrev =
         ((n_par.grid_y / H) .^ tax_prog_scale .* LC) -
         τlev .* ((n_par.grid_y / H) .^ tax_prog_scale .* LC) .^ (1.0 - τprog)
     taxrev[end] =
         n_par.grid_y[end] .* profits -
         τlev .* (n_par.grid_y[end] .* profits) .^ (1.0 - τprog)
-    incgrossaux = ((n_par.grid_y / H) .^ tax_prog_scale .* LC)
+    incgrossaux      = ((n_par.grid_y / H) .^ tax_prog_scale .* LC)
     incgrossaux[end] = n_par.grid_y[end] .* profits
-    av_tax_rate_up = dot(distr_y, taxrev) ./ (dot(distr_y, incgrossaux))
-
-
+    av_tax_rate_up   = dot(distr_y, taxrev) ./ (dot(distr_y, incgrossaux))
 
     # Calculate optimal policies
     # expected margginal values
@@ -279,13 +275,12 @@ function Fsys(
         distrSummaries(PDF_joint, q, c_a_star, c_n_star, n_par, inc, incgross, m_par)
 
     # Error Term on prices/aggregate summary vars (logarithmic, controls)
-    F[indexes.τlev] = av_tax_rate - av_tax_rate_up
-    F[indexes.T] = log(T) - log(dot(distr_y, taxrev) + av_tax_rate * (unionprofits))
-
-    F[indexes.K] = log.(K) - log.(KP)
-    F[indexes.B] = log.(B) - log.(BP)
+    F[indexes.K]  = log.(K) - log.(KP)
+    F[indexes.B]  = log.(B) - log.(BP)
     F[indexes.BD] = log.(BD) - log.(BDact)
     F[indexes.Ht] = log.(Ht) - log.(Htact)
+    F[indexes.τlev] = av_tax_rate - av_tax_rate_up
+    F[indexes.T]  = log(T) - log(dot(distr_y, taxrev) + av_tax_rate * (unionprofits))
 
     # Error Terms on  distribution summaries
     F[indexes.GiniW] = log.(GiniW) - log.(GiniWT)