diff --git a/src/smc_main.jl b/src/smc_main.jl
index 9bc0f96..e8d444d 100644
--- a/src/smc_main.jl
+++ b/src/smc_main.jl
@@ -207,6 +207,20 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
     fixed_para_inds = findall([ θ.fixed for θ in parameters])
     free_para_inds  = findall([!θ.fixed for θ in parameters])
     para_symbols    = [θ.key for θ in parameters]
+    if regime_switching
+        # Concatenate regime symbols for each extra regimes
+        reg_switch_symbols = Vector{Symbol}(undef, n_para_rs)
+        ind = 0
+        for θ in parameters
+            if !isempty(θ.regimes)
+                for i in 2:length(θ.regimes[:value])
+                    ind += 1
+                    reg_switch_symbols[ind] = Symbol(θ.key, "_reg$(i)")
+                end
+            end
+        end
+        push!(para_symbols, reg_switch_symbols...)
+    end
 
     n_free_para     = length(free_para_inds) + n_para_rs
     free_para_inds = vcat(free_para_inds, collect(n_para+1:n_para+n_para_rs))
diff --git a/src/util.jl b/src/util.jl
index 7db4721..9bb191d 100644
--- a/src/util.jl
+++ b/src/util.jl
@@ -134,10 +134,8 @@ function init_stage_print(cloud::Cloud, para_symbols::Vector{Symbol};
     if VERBOSITY[verbose] >= VERBOSITY[:high]
         μ = weighted_mean(cloud)
         σ = weighted_std(cloud)
-        # TODO: need to generalize this output for when parameters are regime-switching. Currently just prints the value
-        # store in p.value
         println("Mean and standard deviation of parameter estimates")
-        for n=1:length(para_symbols)
+        for n = 1:length(para_symbols)
             println("$(para_symbols[n]) = $(round(μ[n], digits = 5)), $(round(σ[n], digits = 5))")
 	    end
     end
@@ -174,8 +172,6 @@ function end_stage_print(cloud::Cloud, para_symbols::Vector{Symbol};
     if VERBOSITY[verbose] >= VERBOSITY[:high]
         μ = weighted_mean(cloud)
         σ = weighted_std(cloud)
-        # TODO: need to generalize this output for when parameters are regime-switching. Currently just prints the value
-        # store in p.value
         println("Mean and standard deviation of parameter estimates")
         for n=1:length(para_symbols)
             println("$(para_symbols[n]) = $(round(μ[n], digits = 5)), $(round(σ[n], digits = 5))")