text:
using FirstOrderOptimization
import StatsBase: sample, WeightVec

# Currently, the gradient descent algorithm seems to converge to the optimum,
# and the code seems to have all the invariances etc expected.
# However, it seems like (for small problem sizes) using the lambda
# from our paper, the true Theta used to generate the data
# does *not* yield the minimal objective value:
# the ThetaHat computed has a *much* lower objective value
# And for large problem sizes even *computing* the objective value (and gradient) take a long time.
# Also, I find it odd that the initialization from the dropping convexity paper
# works as well as it does. Perhaps the gradient computation is wrong?
# Or is 100 x 100 just too small to see any interesting behavior?

# With the nucler norm constrained- (rather than regularized-) problem,
# the code below searches over a range of constraints. It seems like
# the code converges nicely (if slowly).
# When the constraint is set at the nuclear norm of Thetastar,
# the negloglik converges to negloglik(Thetastar).
# But ||ThetaHat - Thetastar||/||Thetastar|| is still more than .1.

################ Now run it

### Generate data

# make it deterministic
srand(10)

m, n, k, ktrue = 300, 300, 2, 1
Theta = randn(m, ktrue) * randn(ktrue, n)
Theta = Theta .- (Theta * ones(n,1)/n)
Theta = Theta / var(Theta)

nobs = 3*m*n #int(10*m*n) # number of observations
K = 4 #int(.05*n) # size of sets
sets = Array(Array{Int,1}, nobs) # set presented
rows = Array(Int, nobs)      # customer who arrived
cols = Array(Int, nobs)      # product they picked
for iobs = 1:nobs
	it = sample(1:m)
	St = sample(1:n, K, replace=false, ordered=true)
	wv = WeightVec(Float64[exp(-Theta[it, j]) for j in St])
	jt = St[sample(wv)]
	rows[iobs], cols[iobs], sets[iobs] = it, jt, St
end
data = MNLdata(nobs, m, n, sets, rows, cols)
println("nucnorm(Theta) = ", nucnorm(Theta))
println("negloglik(Theta) = ", negloglik(Theta, data))

relrmse(xhat, x) = vecnorm(xhat - x) / vecnorm(x) # sqrt(mean((xhat - x).^2)) / sqrt(mean(x.^2))

# per the theorem, we want
# lambda >= 32 sqrt(K*d*log(d)/m/n/nobs)
d = (m+n)/2
lambda = 0

# certify optimality of Theta when solving
# minimize(loss(Theta) + lambda*nucnorm(Theta))
# G is the gradient of the loss at Theta
# lambda is the coefficient in front of the nuclear norm regularization
# here's the idea: (taking lambda = 1 in the discussion, and hoping all nonzero eigenvalues are distinct)
# writing down the optimality conditions of the problem, we see
# 1) on the span (U,V) of the solution Theta,
#    the gradient had better be diagonal with entries all 1
# 2) orthogonal to that span, all we know is that it has to have operator norm <=1
function certify_optimality(G, Theta, lambda; TOL=1e-4)
	# TOL is numerical zero, for our purposes
	u,s,v = svd(Theta)
	ug,sg,vg = svd(G)
	k = sum(s.>=TOL)
	U = u[:,1:k]
	V = v[:,1:k]
	GVVT = (G*V)*V'
	W = (G - U*(U'*G) - GVVT + U*(U'*GVVT))
	notdiag = vecnorm(T + lambda*eye(k)) / k
	notsmall = max(norm(W)-lambda, 0)
	optcrit = notdiag + notsmall
	if optcrit <= TOL
		println("Recovered value is optimal")
	else
		println("Recovered value is not optimal: violates optimality condition
			on the range of the solution by $notdiag
			and off that range by $notsmall")
	end
end

stop(p::SimpleParams, objval, oldobjval, args...; kwargs...) = false

## Frank Wolfe
if false
	println("Fitting a rank $k MNL model with $m rows and $n columns to $nobs observations with choice sets of size $K
	using Frank Wolfe")

	t = zeros(m,n)
	ThetaHat = copy(t)
	@printf("%12s%12s%12s%12s\n", "obj", "delta", "constr", "rel rmse")
	nt = nucnorm(Theta)
	for delta=nt:round(Int, nt/10):nt*2
		# saturate constraint
		if delta > 0 && nucnorm(ThetaHat) > .01
			t = copy(ThetaHat)
			t *= delta/nucnorm(ThetaHat)
		end
		ThetaHat = frank_wolfe(t, delta,
				x->negloglik(x, data),
				x->grad_negloglik(x, data),
				nucnorm, min_lin_st_nucnorm,
				FrankWolfeParams(),
				verbose=false)

		ThetaHatZeroed = ThetaHat .- (ThetaHat * ones(n,1)/n)
		rmset = relrmse(ThetaHatZeroed, Theta)
		@printf("%12.2e%12.2e%12.2e%12.2e\n", negloglik(ThetaHat, data), delta, nucnorm(ThetaHat), rmset)
	end
end

## Prox grad on Theta
if false
	println("Fitting a rank $k MNL model with $m rows and $n columns to $nobs observations with choice sets of size $K
	using the proximal gradient method")
	println("nucnorm(Theta) = ", nucnorm(Theta))
	println("negloglik(Theta) = ", negloglik(Theta, data))

	@show gamma_u = 256*sqrt(K*k*(m+n)*log(m+n)/(m*n*nobs))

	ThetaHat = zeros(m,n)
	@printf("%12s%12s%12s%12s%12s\n", "lambda", "loss", "reg", "obj", "rel rmse")
	nt = nucnorm(Theta)
	for lambda=logspace(1,-4,6)
		ThetaHat = prox_grad!(ThetaHat,
				x->negloglik(x, data),
				x->grad_negloglik(x, data),
				x->lambda*nucnorm(x),
				(x, alpha)->prox_nucnorm(x, alpha*lambda),
				ProxGradParams(ConstantStepSize(gamma_u), 100),
				verbose=true)
		# consider also HopefulStepSize(gamma_u)

		ThetaHatZeroed = ThetaHat .- (ThetaHat * ones(n,1)/n)
		rmset = relrmse(ThetaHatZeroed, Theta)
		loss, reg = negloglik(ThetaHat, data), lambda*nucnorm(ThetaHat)
		@printf("%12.2e%12.2e%12.2e%12.2e%12.2e\n", lambda, loss, reg, loss+reg, rmset)
	end
end

## Factored prox grad
if false
	println("Fitting a rank $k MNL model with $m rows and $n columns to $nobs observations with choice sets of size $K
	using factored proximal gradient descent and varying the regularization")
	println("nucnorm(Theta) = ", nucnorm(Theta))
	println("negloglik(Theta) = ", negloglik(Theta, data))
	@printf("%12s%12s%12s%12s%12s\n", "lambda", "loss", "reg", "obj", "rel rmse")

	# initialize
	lambda0 = 32 * sqrt(K*d*log(d)/m/n/nobs)
	objective = U -> (negloglik(U, data) + lambda0*nucnorm(U))
	grad_objective = U -> (grad_negloglik(U, data) + lambda0 * grad_nucnorm(U))
	U = initialize_dropping_convexity(grad_objective, zeros(m,n), k)

	@show gamma_u = 256*sqrt(K*k*(m+n)*log(m+n)/(m*n*nobs))

	@printf("%12s%12s%12s%12s%12s\n", "lambda", "loss", "reg", "obj", "rel rmse")
	for lambda=logspace(1,-4,6)
		U = prox_grad!(U,
				x->negloglik(x, data),
				x->grad_negloglik(x, data),
				x->lambda*nucnorm(x),
				(x, alpha)->prox_nucnorm(x, alpha*lambda),
				ProxGradParams(HopefulStepSize(), 100),
				verbose=true)
		ThetaHat = rectangular_part(U)
		ThetaHatZeroed = ThetaHat .- (ThetaHat * ones(n,1)/n)
		rmset = relrmse(ThetaHatZeroed, Theta)
		loss, reg = negloglik(ThetaHat, data), lambda*nucnorm(ThetaHat)
		@printf("%12.2e%12.2e%12.2e%12.2e%12.2e\n", lambda, loss, reg, loss+reg, rmset)
	end
end

## Factored gradient descent
if true
	println("Fitting a rank $k MNL model with $m rows and $n columns to $nobs observations with choice sets of size $K
	using factored gradient descent and varying the regularization")
	@printf("%12s%12s%12s%12s%12s\n", "lambda", "loss", "reg", "obj", "rel rmse")

	# initialize
	lambda0 = 32 * sqrt(K*d*log(d)/m/n/nobs)
	objective = U -> (negloglik(U, data) + lambda0*nucnorm(U))
	grad_objective = U -> (grad_negloglik(U, data) + lambda0 * grad_nucnorm(U))
	Uinit = initialize_dropping_convexity(grad_objective, zeros(m,n), k)
	U = copy(Uinit)

	for lambda = lambda0* 2.0^(-10) * 2.0.^(0:0.5:15) # lambda = lambda0*logspace(0,-5,6)
		objective = U -> (negloglik(U, data) + lambda*nucnorm(U))
		grad_objective = U -> (grad_negloglik(U, data) + lambda * grad_nucnorm(U))
		# don't initialize at 0
		if vecnorm(U) < .01*sqrt(k*(m+n))
			U = copy(Uinit)
		end
		U = gradient_descent!(U, objective, grad_objective,
			              SimpleParams(HopefulStepSize(1, .8, 1.5, .1), 200), verbose=false)
		ThetaHat = rectangular_part(U)
		certify_optimality(grad_negloglik(Theta, data), Theta, lambda)
		ThetaHatZeroed = ThetaHat .- (ThetaHat * ones(n,1)/n)
		loss, reg = negloglik(ThetaHat, data), lambda*nucnorm(ThetaHat)
		rmset = relrmse(ThetaHatZeroed, Theta)
		@printf("%12.2e%12.2e%12.2e%12.2e%12.2e\n", lambda, loss, reg, loss+reg, rmset)
	end
end

## gradient descent on Theta

# initialize
# ThetaHat = zeros(m,n)
# ThetaHat = initialize_dropping_convexity(grad_objective, zeros(m,n))

# gradient_descent!(ThetaHat, objective, grad_objective,
# 	                 SimpleParams(DecreasingStepSize(.1), 20), verbose=true)

# compare to the bound if the constant isn't too honkingly large

# we've dropped the constraint Theta*1 = 0.
# the likelihood term is invariant to constants added to each row,
# and the regularizer is minimized when these offsets are chosen to be zero.
# so the solution should approximately satisfy
# println("offsets:\n", ThetaHat*ones(n,1)/n) == 0

# see if we converge to the true solution using
# t = initialize_dropping_convexity(grad_objective, zeros(m,n), k)
# gradient_descent!(t, t->negloglik(t, data), t->grad_negloglik(t, data), SimpleParams(HopefulStepSize(), 50), verbose=true)
# println("relative root mean square error is ", relrmse(rectangular_part(t), Theta))

max_stars_repo_name: madeleineudell/FirstOrderOptimization


text:
"""
`init_subspec!(m::Model1002)`

Initializes a model subspecification by overwriting parameters from
the original model object with new parameter objects. This function is
called from within the model constructor.
"""
function init_subspec!(m::Model1002)
    if subspec(m) == "ss2"
        return
    elseif subspec(m) == "ss8"
        return ss8!(m)
    elseif subspec(m) == "ss9"
        return ss9!(m)
    elseif subspec(m) == "ss10"
        return ss10!(m)
    elseif subspec(m) == "ss11"
        # Normalize s.t. sum of variances of anticipated shocks equals the
        # variance of the contemporaneous shock
        # See measurement equation
        return ss10!(m)
    elseif subspec(m) == "ss12"
        return ss12!(m)
    elseif subspec(m) == "ss13"
        return ss13!(m)
    elseif subspec(m) == "ss14"
        return ss14!(m)
    elseif subspec(m) == "ss15"
        return ss15!(m)
    elseif subspec(m) == "ss16"
        return ss16!(m)
    elseif subspec(m) == "ss17"
        return ss17!(m)
    elseif subspec(m) == "ss18"
        return ss18!(m)
    elseif subspec(m) == "ss19"
        return ss19!(m)
    elseif subspec(m) == "ss20"
        return ss20!(m)
    elseif subspec(m) == "ss51v"
        return ss51v!(m)
    elseif subspec(m) == "ss59"
        return ss59!(m)
    elseif subspec(m) == "ss60"
        return ss60!(m)
    elseif subspec(m) == "ss61"
        return ss61!(m)
    else
        error("This subspec is not defined.")
    end
end

"""
```
ss3!(m::Model1002)
```

"""
function ss3!(m::Model1002)
    ss10!(m)
       m <= parameter(:α, 0.1596, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Normal(0.30, 0.15), fixed=false,
                   description="α: Capital elasticity in the intermediate goods sector's production function (also known as the capital share).",
                   tex_label="\\alpha")

    m <= parameter(:ζ_p, 0.8940, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.3), fixed=false,
                   description="ζ_p: The Calvo parameter. In every period, intermediate goods producers optimize prices with probability (1-ζ_p). With probability ζ_p, prices are adjusted according to a weighted average of the previous period's inflation (π_t1) and steady-state inflation (π_star).",
                   tex_label="\\zeta_p")

    m <= parameter(:ι_p, 0.1865, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.45), fixed=false,
                   description="ι_p: The weight attributed to last period's inflation in price indexation. (1-ι_p) is the weight attributed to steady-state inflation.",
                   tex_label="\\iota_p")

    m <= parameter(:δ, 0.025, fixed=true,
                   description="δ: The capital depreciation rate.",
                   tex_label="\\delta" )

    m <= parameter(:Upsilon, 1.000, (0., 10.), (1e-5, 0.), ModelConstructors.Exponential(), GammaAlt(1., 1.5), fixed=true,
                   description="Υ: The trend evolution of the price of investment goods relative to consumption goods. Set equal to 1.",
                   tex_label="\\Upsilon")

    m <= parameter(:Φ, 1.1066, (1., 10.), (1.00, 10.00), ModelConstructors.Exponential(), Normal(1.25, 0.36), fixed=false,
                   description="Φ: Fixed costs.",
                   tex_label="\\Phi_p")

    m <= parameter(:S′′, 2.7314, (-15., 15.), (-15., 15.), ModelConstructors.Untransformed(), Normal(4., 4.5), fixed=false,
                   description="S'': The second derivative of households' cost of adjusting investment.",
                   tex_label="S''")

    m <= parameter(:h, 0.5347, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Uniform(0., 1.), fixed=false,
                   description="h: Consumption habit persistence.",
                   tex_label="h")

    m <= parameter(:ppsi, 0.6862, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Uniform(0., 1.), fixed=false,
                   description="ppsi: Utilization costs.",
                   tex_label="\\psi")

    m <= parameter(:ν_l, 2.5975, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), Normal(2, 2.25), fixed=false,
                   description="ν_l: The coefficient of relative risk aversion on the labor term of households' utility function.",
                   tex_label="\\nu_l")

    m <= parameter(:ζ_w, 0.9291, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Uniform(0.,1.), fixed=false,
                   description="ζ_w: (1-ζ_w) is the probability with which households can freely choose wages in each period. With probability ζ_w, wages increase at a geometrically weighted average of the steady state rate of wage increases and last period's productivity times last period's inflation.",
                   tex_label="\\zeta_w")

    m <= parameter(:ι_w, 0.2992, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Uniform(0.,1.), fixed=false,
                   description="ι_w: The weight attributed to last period's wage in wage indexation. (1-ι_w) is the weight attributed to steady-state wages.",
                   tex_label="\\iota_w")

    m <= parameter(:λ_w, 1.5000, fixed=true,
                   description="λ_w: The wage markup, which affects the elasticity of substitution between differentiated labor services.",
                   tex_label="\\lambda_w")

    m <= parameter(:β, 0.1402, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), GammaAlt(0.25, 0.3), fixed=false,
                   scaling = x -> 1/(1 + x/100),
                   description="β: Discount rate.",
                   tex_label="100(\\beta^{-1} - 1)")

    m <= parameter(:ψ1, 1.3679, (1e-5, 10.), (1e-5, 10.00), ModelConstructors.Exponential(), Normal(1.5, 0.75), fixed=false,
                   description="ψ₁: Weight on inflation gap in monetary policy rule.",
                   tex_label="\\psi_1")

    m <= parameter(:ψ2, 0.0388, (-0.5, 0.5), (-0.5, 0.5), ModelConstructors.Untransformed(), Normal(0.12, 0.15), fixed=false,
                   description="ψ₂: Weight on output gap in monetary policy rule.",
                   tex_label="\\psi_2")

    m <= parameter(:ψ3, 0.2464, (-0.5, 0.5), (-0.5, 0.5), ModelConstructors.Untransformed(), Normal(0.12, 0.15), fixed=false,
                   description="ψ₃: Weight on rate of change of output gap in the monetary policy rule.",
                   tex_label="\\psi_3")

    m <= parameter(:π_star, 0.5000, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), GammaAlt(0.75, 0.4), fixed=true,
                   scaling = x -> 1 + x/100,
                   description="π_star: The steady-state rate of inflation.",
                   tex_label="\\pi_*")

    m <= parameter(:σ_c, 0.8719, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), Normal(1.5, 1.11), fixed=false,
                   description="σ_c: Coefficient of relative risk aversion.",
                   tex_label="\\sigma_{c}")

    m <= parameter(:ρ, 0.7126, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Uniform(0., 1,), fixed=false,
                   description="ρ: The degree of inertia in the monetary policy rule.",
                   tex_label="\\rho_R")

    m <= parameter(:ϵ_p, 10.000, fixed=true,
                   description="ϵ_p: Curvature parameter in the Kimball aggregator for prices.",
                   tex_label="\\epsilon_{p}")

    m <= parameter(:ϵ_w, 10.000, fixed=true,
                   description="ϵ_w: Curvature parameter in the Kimball aggregator for wages.",
                   tex_label="\\epsilon_{w}")

    # financial frictions parameters
    m <= parameter(:Fω, 0.0300, (1e-5, 0.99999), (1e-5, 0.99), ModelConstructors.SquareRoot(), BetaAlt(0.03, 0.01), fixed=true,
                   scaling = x -> 1 - (1-x)^0.25,
                   description="F(ω): The cumulative distribution function of ω (idiosyncratic iid shock that increases or decreases entrepreneurs' capital).",
                   tex_label="F(\\bar{\\omega})")

    m <= parameter(:spr, 1.7444, (0., 100.), (1e-5, 0.), ModelConstructors.Exponential(), GammaAlt(2., 0.1), fixed=false,
                   scaling = x -> (1 + x/100)^0.25,
                   description="spr_*: Steady-state level of spread.",
                   tex_label="SP_*")

    m <= parameter(:ζ_spb, 0.0559, (1e-5, 0.99999), (1e-5, 0.99), ModelConstructors.SquareRoot(), BetaAlt(0.05, 0.005), fixed=false,
                   description="ζ_spb: The elasticity of the expected exess return on capital (or 'spread') with respect to leverage.",
                   tex_label="\\zeta_{sp,b}")

    m <= parameter(:γ_star, 0.9900, (1e-5, 0.99999), (1e-5, 0.99), ModelConstructors.SquareRoot(), BetaAlt(0.99, 0.002), fixed=true,
                   description="γ_star: Fraction of entrepreneurs who survive and continue operating for another period.",
                   tex_label="\\gamma_*")

    # exogenous processes - level
    m <= parameter(:γ, 0.3673, (-5.0, 5.0), (-5., 5.), ModelConstructors.Untransformed(), Normal(0.4, 0.1), fixed=false,
                   scaling = x -> x/100,
                   description="γ: The log of the steady-state growth rate of technology.",
                   tex_label="100\\gamma")

    m <= parameter(:Lmean, -45.9364, (-1000., 1000.), (-1e3, 1e3), ModelConstructors.Untransformed(), Normal(-45., 5.), fixed=false,
                   description="Lmean: Mean level of hours.",
                   tex_label="\\bar{L}")

    m <= parameter(:g_star, 0.1800, fixed=true,
                   description="g_star: 1 - (c_star + i_star)/y_star.",
                   tex_label="g_*")

    # exogenous processes - autocorrelation
    m <= parameter(:ρ_g, 0.9863, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Uniform(0., 1.), fixed=false,
                   description="ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label="\\rho_g")

    m <= parameter(:ρ_b, 0.9410, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Uniform(0., 1.), fixed=false,
                   description="ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
                   tex_label="\\rho_b")

    m <= parameter(:ρ_μ, 0.8735, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Uniform(0., 1.), fixed=false,
                   description="ρ_μ: AR(1) coefficient in capital adjustment cost process.",
                   tex_label="\\rho_{\\mu}")

    m <= parameter(:ρ_ztil, 0.9446, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Uniform(0., 1.), fixed=false,
                   description="ρ_ztil: AR(1) coefficient in the technology process.",
                   tex_label="\\rho_{\\tilde{z}}")

    m <= parameter(:ρ_λ_f, 0.8827, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Uniform(0., 1.), fixed=false,
                   description="ρ_λ_f: AR(1) coefficient in the price mark-up shock process.",
                   tex_label="\\rho_{\\lambda_f}")

    m <= parameter(:ρ_λ_w, 0.3884, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Uniform(0., 1.), fixed=false,
                   description="ρ_λ_w: AR(1) coefficient in the wage mark-up shock process.",
                   tex_label="\\rho_{\\lambda_w}")

    # monetary policy shock - see eqcond
    m <= parameter(:ρ_rm, 0.2135, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Uniform(0., 1.), fixed=false,
                   description="ρ_rm: AR(1) coefficient in the monetary policy shock process.",
                   tex_label="\\rho_{r^m}")

    m <= parameter(:ρ_σ_w, 0.9898, (1e-5, 0.99999), (1e-5, 0.99), ModelConstructors.SquareRoot(), Uniform(0., 1.0), fixed=false,
                   description="ρ_σ_w: The standard deviation of entrepreneurs' capital productivity follows an exogenous process with mean ρ_σ_w. Innovations to the process are called _spread shocks_.",
                   tex_label="\\rho_{\\sigma_\\omega}")

    m <= parameter(:ρ_μ_e, 0.7500, (1e-5, 0.99999), (1e-5, 0.99), ModelConstructors.SquareRoot(), Uniform(0., 1.0), fixed=true,
                   description="ρ_μ_e: AR(1) coefficient in the exogenous bankruptcy cost process.",
                   tex_label="\\rho_{\\mu_e}")

    m <= parameter(:ρ_γ, 0.7500, (1e-5, 0.99999), (1e-5, 0.99), ModelConstructors.SquareRoot(), Uniform(0., 1.0), fixed=true,
                   description="ρ_γ: AR(1) coefficient in the process describing the fraction of entrepreneurs surviving period t.",
                   tex_label="\\rho_{\\gamma}")

    m <= parameter(:ρ_π_star, 0.9900, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Uniform(0., 1.), fixed=true,
                   description="ρ_π_star: AR(1) coefficient in the time-varying inflation target process.",
                   tex_label="\\rho_{\\pi_*}")

    m <= parameter(:ρ_lr, 0.6936, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Uniform(0., 1.), fixed=false,
                   tex_label="\\rho_{10y}")

    m <= parameter(:ρ_z_p, 0.8910, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Uniform(0., 1.), fixed=false,
                   description="ρ_z_p: AR(1) coefficient in the process describing the permanent component of productivity.",
                   tex_label="\\rho_{z^p}")

    m <= parameter(:ρ_tfp, 0.1953, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Uniform(0., 1.), fixed=false,
                   tex_label="\\rho_{tfp}")

    m <= parameter(:ρ_gdpdef, 0.5379, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Uniform(0., 1.), fixed=false,
                   tex_label="\\rho_{gdpdef}")

    m <= parameter(:ρ_corepce, 0.2320, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Uniform(0., 1.), fixed=false,
                   tex_label="\\rho_{pce}")

    m <= parameter(:ρ_gdp, 0., (-0.999, 0.999), (-0.999, 0.999), ModelConstructors.SquareRoot(), Normal(0.0, 0.2), fixed=false,
                   tex_label="\\rho_{gdp}")

    m <= parameter(:ρ_gdi, 0., (-0.999, 0.999), (-0.999, 0.999), ModelConstructors.SquareRoot(), Normal(0.0, 0.2), fixed=false,
                   tex_label="\\rho_{gdi}")

    m <= parameter(:ρ_gdpvar, 0., (-0.999, 0.999), (-0.999, 0.999), ModelConstructors.SquareRoot(), Normal(0.0, 0.4), fixed=false,
                   tex_label="\\varrho_{gdp}")

    m <= parameter(:me_level, 1., (-0.999, 0.999), (-0.999, 0.999), ModelConstructors.Untransformed(), Normal(0.0, 0.4), fixed=true,
                   description="me_level: Indicator of cointegration of GDP and GDI.",
                   tex_label="\\mathcal{C}_{me}")

    # exogenous processes - standard deviation
    m <= parameter(:σ_g, 2.5230, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_g: The standard deviation of the government spending process.",
                   tex_label="\\sigma_{g}")

    m <= parameter(:σ_b, 0.0292, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_b: The standard deviation of the intertemporal preference shifter process.",
                   tex_label="\\sigma_{b}")

    m <= parameter(:σ_μ, 0.4559, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
                   tex_label="\\sigma_{\\mu}")

    m <= parameter(:σ_ztil, 0.6742, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_ztil: The standard deviation of the process describing the stationary component of productivity.",
                   tex_label="\\sigma_{\\tilde{z}}")

    m <= parameter(:σ_λ_f, 0.1314, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_λ_f: The mean of the process that generates the price elasticity of the composite good. Specifically, the elasticity is (1+λ_{f,t})/(λ_{f_t}).",
                   tex_label="\\sigma_{\\lambda_f}")

    m <= parameter(:σ_λ_w, 0.3864, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{\\lambda_w}")

    m <= parameter(:σ_r_m, 0.2380, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_r_m: The standard deviation of the monetary policy shock.",
                   tex_label="\\sigma_{r^m}")

    m <= parameter(:σ_σ_ω, 0.0428, (1e-7,100.), (1e-5, 0.), ModelConstructors.Exponential(), RootInverseGamma(4, 0.05), fixed=false,
                   description="σ_σ_ω: The standard deviation of entrepreneurs' capital productivity follows an exogenous process with standard deviation σ_σ_ω.",
                   tex_label="\\sigma_{\\sigma_\\omega}")

    m <= parameter(:σ_μ_e, 0.0000, (0., 100.), (1e-5, 0.), ModelConstructors.Exponential(), RootInverseGamma(4, 0.05), fixed=true,
                   description="σ_μ_e: Exogenous bankrupcy costs follow an exogenous process with standard deviation σ_μ_e.",
                   tex_label="\\sigma_{\\mu_e}")

    m <= parameter(:σ_γ, 0.0000, (0., 100.), (1e-5, 0.), ModelConstructors.Exponential(), RootInverseGamma(4, 0.01), fixed=true,
                   description="σ_γ: The fraction of entrepreneurs surviving period t follows an exogenous process with standard deviation σ_γ.",
                   tex_label="\\sigma_{\\gamma}")

    m <= parameter(:σ_π_star, 0.0269, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(6, 0.03), fixed=false,
                   description="σ_π_star: The standard deviation of the inflation target.",
                   tex_label="\\sigma_{\\pi_*}")

    m <= parameter(:σ_lr, 0.1766, (1e-8,10.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.75), fixed=false,
                   tex_label="\\sigma_{10y}")

    m <= parameter(:σ_z_p, 0.1662, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_z_p: The standard deviation of the shock to the permanent component of productivity.",
                   tex_label="\\sigma_{z^p}")

    m <= parameter(:σ_tfp, 0.9391, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{tfp}")

    m <= parameter(:σ_gdpdef, 0.1575, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{gdpdef}")

    m <= parameter(:σ_corepce, 0.0999, (1e-8, 5.),(1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{pce}")

    m <= parameter(:σ_gdp, 0.1, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{gdp}")

    m <= parameter(:σ_gdi, 0.1, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{gdi}")

  # standard deviations of the anticipated policy shocks
    for i = 1:n_mon_anticipated_shocks_padding(m)
        if i <= n_mon_anticipated_shocks(m)
            m <= parameter(Symbol("σ_r_m$i"), .2, (1e-7, 100.), (1e-5, 0.), ModelConstructors.Exponential(), RootInverseGamma(4, .2), fixed=false,
                           description="σ_r_m$i: Standard deviation of the $i-period-ahead anticipated policy shock.",
                           tex_label=@sprintf("\\sigma_{ant%d}",i))
        else
            m <= parameter(Symbol("σ_r_m$i"), .0, (1e-7, 100.), (1e-5, 0.), ModelConstructors.Exponential(), RootInverseGamma(4, .2), fixed=true,
                           description="σ_r_m$i: Standard deviation of the $i-period-ahead anticipated policy shock.",
                           tex_label=@sprintf("\\sigma_{ant%d}",i))
        end
    end

    m <= parameter(:η_gz, 0.8400, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.50, 0.20), fixed=false,
                   description="η_gz: Correlate g and z shocks.",
                   tex_label="\\eta_{gz}")

    m <= parameter(:η_λ_f, 0.7892, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.50, 0.20), fixed=false,
                   description="η_λ_f: Moving average component in the price markup shock.",
                   tex_label="\\eta_{\\lambda_f}")

    m <= parameter(:η_λ_w, 0.4226, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.50, 0.20), fixed=false,
                   description="η_λ_w: Moving average component in the wage markup shock.",
                   tex_label="\\eta_{\\lambda_w}")

    m <= parameter(:Iendoα, 0.0000, (0.000, 1.000), (0., 0.), ModelConstructors.Untransformed(), BetaAlt(0.50, 0.20), fixed=true,
                   description="Iendoα: Indicates whether to use the model's endogenous α in the capacity utilization adjustment of total factor productivity.",
                   tex_label="I\\{\\alpha^{model}\\}")

    m <= parameter(:Γ_gdpdef, 1.0354, (-10., 10.), (-10., -10.), ModelConstructors.Untransformed(), Normal(1.00, 2.), fixed=false,
                   tex_label="\\gamma_{gdpdef}")

    m <= parameter(:δ_gdpdef, 0.0181, (-10., 10.), (-10., -10.), ModelConstructors.Untransformed(), Normal(0.00, 2.), fixed=false,
                   tex_label="\\delta_{gdpdef}")

    m <= parameter(:γ_gdi, 1., (-10., 10.), (-10., -10.), ModelConstructors.Untransformed(), Normal(1., 2.), fixed=false,
                   tex_label="\\gamma_{gdi}")

    m <= parameter(:δ_gdi, 0., (-10., 10.), (-10., -10.), ModelConstructors.Untransformed(), Normal(0.00, 2.), fixed=false,
                   tex_label="\\delta_{gdi}")

end

"""
```
ss5!(m::Model1002)
```
The loose prior subspec we ACTUALLY WANT since it's only loose on the parameters we allow to switch
"""
function ss5!(m::Model1002)
    ss10!(m)
    m <= parameter(:ψ1, 1.3679, (1e-5, 10.), (1e-5, 10.00), ModelConstructors.Exponential(), Normal(1.5, 0.75), fixed=false,
                   description="ψ₁: Weight on inflation gap in monetary policy rule.",
                   tex_label="\\psi_1")

    m <= parameter(:ψ2, 0.0388, (-0.5, 0.5), (-0.5, 0.5), ModelConstructors.Untransformed(), Normal(0.12, 0.15), fixed=false,
                   description="ψ₂: Weight on output gap in monetary policy rule.",
                   tex_label="\\psi_2")

    m <= parameter(:ψ3, 0.2464, (-0.5, 0.5), (-0.5, 0.5), ModelConstructors.Untransformed(), Normal(0.12, 0.15), fixed=false,
                   description="ψ₃: Weight on rate of change of output gap in the monetary policy rule.",
                   tex_label="\\psi_3")

    m <= parameter(:ρ, 0.7126, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.75, 0.30), fixed=false,
                   description="ρ: The degree of inertia in the monetary policy rule.",
                   tex_label="\\rho_R")

    m <= parameter(:ζ_p, 0.8940, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.3), fixed=false,
                   description="ζ_p: The Calvo parameter (Regime 1). In every period, intermediate goods producers optimize prices with probability (1-ζ_p). With probability ζ_p, prices are adjusted according to a weighted average of the previous period's inflation (π_t1) and steady-state inflation (π_star).",
                   tex_label="\\zeta_p")
end

"""
```
ss8!(m::Model1002)
```

Initializes subspec 8 of `Model1002`. In this subspecification, γ_gdi
and δ_gdi are fixed at 1 and 0 respectively.
"""
function ss8!(m::Model1002)
    m <= parameter(:γ_gdi, 1., (-10., 10.), (-10., -10.), Untransformed(), Normal(1., 2.), fixed = true,
                   tex_label = "\\gamma_{gdi}")

    m <= parameter(:δ_gdi, 0., (-10., 10.), (-10., -10.), Untransformed(), Normal(0.00, 2.), fixed = true,
                   tex_label = "\\delta_{gdi}")
end

"""
```
ss9!(m::Model1002)
```

Initializes subspec 9 of `Model1002`. This subspecification is ss8 + Iskander's
changes + the bounds for beta-distributed parameters have been changed from
(1e-5, 0.99) to (0.0, 1.0).
"""
function ss9!(m::Model1002)

    ss8!(m)

    m <= parameter(:ζ_p, 0.8940, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.1), fixed = false,
                   description = "ζ_p: The Calvo parameter. In every period, intermediate goods producers optimize prices with probability (1-ζ_p). With probability ζ_p, prices are adjusted according to a weighted average of the previous period's inflation (π_t1) and steady-state inflation (π_star).",
                   tex_label = "\\zeta_p")

    m <= parameter(:ι_p, 0.1865, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.15), fixed = false,
                   description = "ι_p: The weight attributed to last period's inflation in price indexation. (1-ι_p) is the weight attributed to steady-state inflation.",
                   tex_label = "\\iota_p")

    m <= parameter(:h, 0.5347, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.7, 0.1), fixed = false,
                   description = "h: Consumption habit persistence.",
                   tex_label = "h")

    m <= parameter(:ppsi, 0.6862, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.15), fixed = false,
                   description = "ppsi: Utilization costs.",
                   tex_label = "\\psi")

    m <= parameter(:ζ_w, 0.9291, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.1), fixed = false,
                   description = "ζ_w: (1-ζ_w) is the probability with which households can freely choose wages in each period. With probability ζ_w, wages increase at a geometrically weighted average of the steady state rate of wage increases and last period's productivity times last period's inflation.",
                   tex_label = "\\zeta_w")

    m <= parameter(:ι_w, 0.2992, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.15), fixed = false,
                   description = "ι_w: The weight attributed to last period's wage in wage indexation. (1-ι_w) is the weight attributed to steady-state wages.",
                   tex_label = "\\iota_w")

    m <= parameter(:ρ, 0.7126, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.75, 0.10), fixed = false,
                   description = "ρ: The degree of inertia in the monetary policy rule.",
                   tex_label = "\\rho_R")

    # Financial frictions parameters
    m <= parameter(:Fω, 0.0300, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.03, 0.01), fixed = true,
                   scaling = x -> 1 - (1-x)^0.25,
                   description = "F(ω): The cumulative distribution function of ω (idiosyncratic iid shock that increases or decreases entrepreneurs' capital).",
                   tex_label = "F(\\bar{\\omega})")

    m <= parameter(:ζ_spb, 0.0559, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.05, 0.005), fixed = false,
                   description = "ζ_spb: The elasticity of the expected exess return on capital (or 'spread') with respect to leverage.",
                   tex_label = "\\zeta_{sp,b}")

    m <= parameter(:γ_star, 0.9900, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.99, 0.002),
                   fixed = true,
                   description = "γ_star: Fraction of entrepreneurs who survive and continue operating for another period.",
                   tex_label = "\\gamma_*")

    m <= parameter(:ρ_g, 0.9863, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label = "\\rho_g")

    m <= parameter(:ρ_b, 0.9410, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
                   tex_label = "\\rho_b")

    m <= parameter(:ρ_μ, 0.8735, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_μ: AR(1) coefficient in capital adjustment cost process.",
                   tex_label = "\\rho_{\\mu}")

    m <= parameter(:ρ_ztil, 0.9446, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_ztil: AR(1) coefficient in the technology process.",
                   tex_label = "\\rho_ztil")

    m <= parameter(:ρ_λ_f, 0.8827, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_λ_f: AR(1) coefficient in the price mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_f}")

    m <= parameter(:ρ_λ_w, 0.3884, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_λ_w: AR(1) coefficient in the wage mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_w}")

    m <= parameter(:ρ_rm, 0.2135, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_rm: AR(1) coefficient in the monetary policy shock process.",
                   tex_label = "\\rho_{r^m}")

    m <= parameter(:ρ_σ_w, 0.9898, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.75, 0.15), fixed = false,
                   description = "ρ_σ_w: The standard deviation of entrepreneurs' capital productivity follows an exogenous process with mean ρ_σ_w. Innovations to the process are called _spread shocks_.",
                   tex_label = "\\rho_{\\sigma_\\omega}")

    m <= parameter(:ρ_μ_e, 0.7500, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.75, 0.15), fixed = true,
                   description = "ρ_μ_e: AR(1) coefficient in the exogenous bankruptcy cost process.",
                   tex_label = "\\rho_{\\mu_e}")

    m <= parameter(:ρ_γ, 0.7500, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.75, 0.15), fixed = true,
                   description = "ρ_γ: AR(1) coefficient in the process describing the fraction of entrepreneurs surviving period t.",
                   tex_label = "\\rho_{\\gamma}")

    m <= parameter(:ρ_π_star, 0.9900, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_π_star: AR(1) coefficient in the time-varying inflation target process.",
                   tex_label = "\\rho_{\\pi_*}")

    m <= parameter(:ρ_lr, 0.6936, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   tex_label = "\\rho_{10y}")

    m <= parameter(:ρ_z_p, 0.8910, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_z_p: AR(1) coefficient in the process describing the permanent component of productivity.",
                   tex_label = "\\rho_{z^p}")

    m <= parameter(:ρ_tfp, 0.1953, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   tex_label = "\\rho_{tfp}")

    m <= parameter(:ρ_gdpdef, 0.5379, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   tex_label = "\\rho_{gdpdef}")

    m <= parameter(:ρ_corepce, 0.2320, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   tex_label = "\\rho_{pce}")

    m <= parameter(:η_gz, 0.8400, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.50, 0.20), fixed = false,
                   description = "η_gz: Correlate g and z shocks.",
                   tex_label = "\\eta_{gz}")

    m <= parameter(:η_λ_f, 0.7892, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.50, 0.20), fixed = false,
                   tex_label = "\\eta_{\\lambda_f}")

    m <= parameter(:η_λ_w, 0.4226, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.50, 0.20), fixed = false,
                   description = "η_λ_w: AR(2) coefficient on wage markup shock process.",
                   tex_label = "\\eta_{\\lambda_w}")

    m <= parameter(:Iendoα, 0.0000, (0.0, 1.0), (0.0, 0.0), Untransformed(), BetaAlt(0.50, 0.20), fixed = true,
                   description = "Iendoα: Indicates whether to use the model's endogenous α in the capacity utilization adjustment of total factor productivity.",
                   tex_label = "I\\{\\alpha^{model}\\}")
end

"""
```
ss10!(m::Model1002)
```

Initializes subspec 10 of `Model1002`. This subspecification is the same as ss9,
except `betabar` is defined with `m[:σ_c]` instead of `σ_ω_star`.
"""
function ss10!(m::Model1002)
    ss9!(m)
end

"""
```
ss12!(m::Model1002)
```

Initializes subspec 12 of `Model1002`. This subspecification is the same as ss10,
but exogenous processes are added as pseudo-observables.
"""
function ss12!(m::Model1002)
    ss10!(m)
end

"""
```
ss13!(m::Model1002)
```

Initializes subspec 13 of `Model1002`. This subspecification is the same as ss10,
but tracks Sinf_t, the contribution of expected future marginal cost, as a pseudo-obs.
"""
function ss13!(m::Model1002)
    ss10!(m)
end

"""
```
ss14!(m::Model1002)
```

Initializes subspec 14 of `Model1002`. This subspecification is the same as ss13,
but with stationarity in the levels (rather than growth) of the measurement
error in the TFP process.
"""
function ss14!(m::Model1002)
    ss13!(m)
end

"""
```
ss15!(m::Model1002)
```

Initializes subspec 15 of `Model1002`. This subspecification is the same as ss14,
but with the TFP observables series now being Fernald's TFP adjusted to utilization.
"""
function ss15!(m::Model1002)
    ss14!(m)
end

"""
```
ss16!(m::Model1002)
```

Initializes subspec 16 of `Model1002`. This subspecification is the same as ss15,
but with the measurement equation using log labor share instead of wage growth.
"""
function ss16!(m::Model1002)
    ss15!(m)
end

"""
```
ss17!(m::Model1002)
```

Initializes subspec 17 of `Model1002`. This subspecification is the same as ss13,
but with the measurement equation using log labor share instead of wage growth.
"""
function ss17!(m::Model1002)
    ss13!(m)
end

"""
```
ss18!(m::Model1002)
```

Initializes subspec 18 of `Model1002`. This subspecification is the same as ss14,
but with ρ_tfp = 0 (fixed) and a near-zero prior on σ_tfp
"""
function ss18!(m::Model1002)
    ss14!(m)
    m <= parameter(:ρ_tfp, 0., (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true, tex_label = "\\rho_{tfp}")
    m <= parameter(:σ_tfp, 0.01, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, sqrt(0.0001)), fixed=false,
                   tex_label="\\sigma_{tfp}")
end

"""
```
ss19!(m::Model1002)
```

Initializes subspec 19 of `Model1002`. This subspecification is the same as ss14,
but with a near zero prior on ρ_tfp and σ_tfp
"""
function ss19!(m::Model1002)
    ss14!(m)
    m <= parameter(:ρ_tfp, 0.01, (1e-5, 0.999999), (1e-5, 0.999999), ModelConstructors.SquareRoot(),
                   BetaAlt(0.01, 0.02), fixed = false, tex_label = "\\rho_{tfp}")
    m <= parameter(:σ_tfp, 0.01, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, sqrt(0.0001)), fixed=false,
                   tex_label="\\sigma_{tfp}")
end

"""
```
ss20!(m::Model1002)
```

Initializes subspec 20 of `Model1002`. This subspecification is the same as ss10,
but the coefficient on markup shocks in the price Phillips curve is re-scaled.
"""
function ss20!(m::Model1002)
    ss10!(m)
end

"""
```
ss51v!(m::Model1002)
```

Second regime where b, mp, and anticipated shocks are ctive.
"""
function ss51v!(m::Model1002)
    ModelConstructors.set_regime_val!(m[:σ_g], 2, 0., override_bounds = true)
    ModelConstructors.set_regime_val!(m[:σ_b], 2, 0.0292)
    ModelConstructors.set_regime_val!(m[:σ_μ], 2, 0., override_bounds = true)
    ModelConstructors.set_regime_val!(m[:σ_ztil], 2, 0., override_bounds = true)
    ModelConstructors.set_regime_val!(m[:σ_λ_f], 2, 0., override_bounds = true)
    ModelConstructors.set_regime_val!(m[:σ_λ_w], 2, 0., override_bounds = true)
    ModelConstructors.set_regime_val!(m[:σ_r_m], 2, 0., override_bounds = true)
    ModelConstructors.set_regime_val!(m[:σ_σ_ω], 2, 0., override_bounds = true)
    ModelConstructors.set_regime_val!(m[:σ_μ_e], 2, 0., override_bounds = true)
    ModelConstructors.set_regime_val!(m[:σ_γ], 2, 0., override_bounds = true)
    ModelConstructors.set_regime_val!(m[:σ_π_star], 2, 0., override_bounds = true)
    ModelConstructors.set_regime_val!(m[:σ_lr], 2, 0., override_bounds = true)
    ModelConstructors.set_regime_val!(m[:σ_z_p], 2, 0., override_bounds = true)
    ModelConstructors.set_regime_val!(m[:σ_tfp], 2, 0., override_bounds = true)
    ModelConstructors.set_regime_val!(m[:σ_gdpdef], 2, 0., override_bounds = true)
    ModelConstructors.set_regime_val!(m[:σ_corepce], 2, 0., override_bounds = true)
    ModelConstructors.set_regime_val!(m[:σ_gdp], 2, 0., override_bounds = true)
    ModelConstructors.set_regime_val!(m[:σ_gdi], 2, 0., override_bounds = true)

    # standard deviations of the anticipated policy shocks
    for i = 1:n_mon_anticipated_shocks(m)
        ModelConstructors.set_regime_val!(m[Symbol("σ_r_m$(i)")], 2, 0., override_bounds = true)
    end

    # Add regime switching settings
    m <= Setting(:regime_switching, true)
    m <= Setting(:n_regimes, 1)
    m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m))) # By default, 1 regime starting at the presample start date.
end

"""
```
ss59!(m::Model1002)
```

models a "temporary shutdown" scenario of the COVID-19 pandemic by using the following features:

1. biidc, φ, and ziid shocks alive in 2020Q1-Q2
2. One-period ahead anticipated shocks for the iid shocks alive in 2020Q1 and expected to hit in 2020Q2.
     These shocks are assumed to be twice the size of the unanticipated shocks which are realized in 2020Q1.
3. The standard deviation of the zp (permanent technology shock) is zero.
4. The standard deviation of the other shocks are set to 0.25 times their value in the time periods before 2020Q1.
5. The standard deviation of monetary policy shocks are set to 10 times their value in the time periods before 2020Q1.

Confer with the official documentation about the COVID-19 shocks on the GitHub page.
"""
function ss59!(m::Model1002)
    # Set up regime-switching
    m <= Setting(:regime_switching, true)
    m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2020, 3, 31),
                                                3 => Date(2020, 6, 30), 4 => Date(2020, 9, 30)))
    setup_regime_switching_inds!(m)

    # Allow the lower bound of non-COVID-19 parameters to equal zero
    m <= parameter(:σ_g, 2.5230, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_g: The standard deviation of the government spending process.",
                   tex_label="\\sigma_{g}")

    m <= parameter(:σ_b, 0.0292, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_b: The standard deviation of the intertemporal preference shifter process.",
                   tex_label="\\sigma_{b}")

    m <= parameter(:σ_μ, 0.4559, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
                   tex_label="\\sigma_{\\mu}")

    m <= parameter(:σ_ztil, 0.6742, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_ztil: The standard deviation of the process describing the stationary component of productivity.",
                   tex_label="\\sigma_{\\tilde{z}}")

    m <= parameter(:σ_λ_f, 0.1314, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_λ_f: The mean of the process that generates the price elasticity of the composite good. Specifically, the elasticity is (1+λ_{f,t})/(λ_{f_t}).",
                   tex_label="\\sigma_{\\lambda_f}")

    m <= parameter(:σ_λ_w, 0.3864, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{\\lambda_w}")

    m <= parameter(:σ_r_m, 0.2380, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_r_m: The standard deviation of the monetary policy shock.",
                   tex_label="\\sigma_{r^m}")

    m <= parameter(:σ_σ_ω, 0.0428, (0., 100.), (1e-5, 0.), ModelConstructors.Exponential(), RootInverseGamma(4, 0.05), fixed=false,
                   description="σ_σ_ω: The standard deviation of entrepreneurs' capital productivity follows an exogenous process with standard deviation σ_σ_ω.",
                   tex_label="\\sigma_{\\sigma_\\omega}")

    m <= parameter(:σ_μ_e, 0.0000, (0. ,100.), (1e-5, 0.), ModelConstructors.Exponential(), RootInverseGamma(4, 0.05), fixed=false,
                   description="σ_μ_e: Exogenous bankrupcy costs follow an exogenous process with standard deviation σ_μ_e.",
                   tex_label="\\sigma_{\\mu_e}")

    m <= parameter(:σ_γ, 0.0000, (0.,100.), (1e-5, 0.), ModelConstructors.Exponential(), RootInverseGamma(4, 0.01), fixed=false,
                   description="σ_γ: The fraction of entrepreneurs surviving period t follows an exogenous process with standard deviation σ_γ.",
                   tex_label="\\sigma_{\\gamma}")

    m <= parameter(:σ_π_star, 0.0269, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(6, 0.03), fixed=false,
                   description="σ_π_star: The standard deviation of the inflation target.",
                   tex_label="\\sigma_{\\pi_*}")

    m <= parameter(:σ_lr, 0.1766, (0.,10.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.75), fixed=false,
                   tex_label="\\sigma_{10y}")

    m <= parameter(:σ_z_p, 0.1662, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_z_p: The standard deviation of the shock to the permanent component of productivity.",
                   tex_label="\\sigma_{z^p}")

    m <= parameter(:σ_tfp, 0.9391, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{tfp}")

    m <= parameter(:σ_gdpdef, 0.1575, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{gdpdef}")

    m <= parameter(:σ_corepce, 0.0999, (0., 5.),(0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{pce}")

    m <= parameter(:σ_gdp, 0.1, (0., 5.),(0., 5.),ModelConstructors.Exponential(),RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{gdp}")

    m <= parameter(:σ_gdi, 0.1, (0., 5.),(0., 5.),ModelConstructors.Exponential(),RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{gdi}")

    # standard deviations of the anticipated policy shocks
    for i = 1:n_mon_anticipated_shocks(m)
        m <= parameter(Symbol("σ_r_m$(i)"), .2, (0., 100.), (1e-5, 0.), ModelConstructors.Exponential(),
                       RootInverseGamma(4, .2), fixed=false,
                       description="σ_r_m$(i): Standard deviation of the $i-period-ahead anticipated policy shock.",
                       tex_label=@sprintf("\\sigma_{ant%d}",i))
    end
    for i = n_mon_anticipated_shocks(m) + 1:n_mon_anticipated_shocks_padding(m)
        m <= parameter(Symbol("σ_r_m$(i)"), 0., (0., 100.), (1e-5, 0.), ModelConstructors.Exponential(),
                       RootInverseGamma(4, .2), fixed=false,
                       description="σ_r_m$(i): Standard deviation of the $i-period-ahead anticipated policy shock.",
                       tex_label=@sprintf("\\sigma_{ant%d}",i))
    end

    # Define regimes for standard shocks
    for i in 1:4
        adj = (i == 2 || i == 3) ? .25 : 1.
        ModelConstructors.set_regime_val!(m[:σ_g], i, adj * m[:σ_g].value)
        ModelConstructors.set_regime_val!(m[:σ_b], i, adj * m[:σ_b].value)
        ModelConstructors.set_regime_val!(m[:σ_μ], i, adj * m[:σ_μ].value)
        ModelConstructors.set_regime_val!(m[:σ_ztil], i, adj * m[:σ_ztil].value)
        ModelConstructors.set_regime_val!(m[:σ_λ_f], i, adj * m[:σ_λ_f].value)
        ModelConstructors.set_regime_val!(m[:σ_λ_w], i, adj * m[:σ_λ_w].value)
        ModelConstructors.set_regime_val!(m[:σ_r_m], i, adj * m[:σ_r_m].value)
        ModelConstructors.set_regime_val!(m[:σ_σ_ω], i, adj * m[:σ_σ_ω].value)
        ModelConstructors.set_regime_val!(m[:σ_μ_e], i, adj * m[:σ_μ_e].value)
        ModelConstructors.set_regime_val!(m[:σ_γ], i, adj * m[:σ_γ].value)
        ModelConstructors.set_regime_val!(m[:σ_π_star], i, adj * m[:σ_π_star].value)
        ModelConstructors.set_regime_val!(m[:σ_lr], i, adj * m[:σ_lr].value)
        ModelConstructors.set_regime_val!(m[:σ_z_p], i, adj * m[:σ_z_p].value)
        ModelConstructors.set_regime_val!(m[:σ_tfp], i, adj * m[:σ_tfp].value)
        err_adj = (i == 3) ? 10. : adj
        ModelConstructors.set_regime_val!(m[:σ_gdpdef], i, err_adj * m[:σ_gdpdef].value)
        ModelConstructors.set_regime_val!(m[:σ_corepce], i, err_adj * m[:σ_corepce].value)
        ModelConstructors.set_regime_val!(m[:σ_gdp], i, adj * m[:σ_gdp].value)
        ModelConstructors.set_regime_val!(m[:σ_gdi], i, adj * m[:σ_gdi].value)

        for j = 1:DSGE.n_mon_anticipated_shocks(m)
            ModelConstructors.set_regime_val!(m[Symbol("σ_r_m$(j)")], i, adj * m[Symbol("σ_r_m$(j)")].value)
        end
    end

    ModelConstructors.set_regime_val!(m[:σ_r_m], 2, 10. * m[:σ_r_m].value) # By default, we set these to 10 times their
    ModelConstructors.set_regime_val!(m[:σ_r_m], 3, 10. * m[:σ_r_m].value) # value in the pre-COVID and post-COVID regimes
    for i = 2:3
        for j = 1:n_mon_anticipated_shocks(m)
            ModelConstructors.set_regime_val!(m[Symbol("σ_r_m$(j)")], i, m[Symbol("σ_r_m$(j)")].value)
        end
    end
    ModelConstructors.set_regime_val!(m[:σ_z_p], 2, 0.)
    ModelConstructors.set_regime_val!(m[:σ_z_p], 3, 0.)

    # Set up the COVID-19 shocks. We do not add anticipation by default
    ModelConstructors.set_regime_val!(m[:σ_φ], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_φ], 2, 400.) # initial value from literature, blown up by 100
    ModelConstructors.set_regime_val!(m[:σ_φ], 3, 400.)
    ModelConstructors.set_regime_val!(m[:σ_φ], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_φ_prop], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_φ_prop], 2, 2.)
    ModelConstructors.set_regime_val!(m[:σ_φ_prop], 3, 0.)
    ModelConstructors.set_regime_val!(m[:σ_φ_prop], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc], 2, 4.)
    ModelConstructors.set_regime_val!(m[:σ_biidc], 3, 4.)
    ModelConstructors.set_regime_val!(m[:σ_biidc], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc_prop], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc_prop], 2, 2.)
    ModelConstructors.set_regime_val!(m[:σ_biidc_prop], 3, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc_prop], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid], 2, 5.)
    ModelConstructors.set_regime_val!(m[:σ_ziid], 3, 5.)
    ModelConstructors.set_regime_val!(m[:σ_ziid], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid_prop], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid_prop], 2, 2.)
    ModelConstructors.set_regime_val!(m[:σ_ziid_prop], 3, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid_prop], 4, 0.)
end

"""
```
ss60!(m::Model1002)
```
models a "temporary shutdown with business cycle dynamics" scenario of the COVID-19 pandemic by using the following features:

1. biidc, φ, and ziid shocks alive in 2020Q1-Q2
2. One-period ahead anticipated shocks for the iid shocks alive in 2020Q1 and expected to hit in 2020Q2.
     These shocks are assumed to be twice the size of the unanticipated shocks which are realized in 2020Q1.
3. The standard deviation of monetary policy shocks are set to 10 times their value in the time periods before 2020Q1.

Confer with the official documentation about the COVID-19 shocks on the GitHub page.
"""
function ss60!(m::Model1002)
    # Set up regime-switching
    m <= Setting(:regime_switching, true)
    m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2020, 3, 31),
                                                3 => Date(2020, 6, 30), 4 => Date(2020, 9, 30)))
    setup_regime_switching_inds!(m)

    # Allow the lower bound of non-COVID-19 parameters to equal zero
    m <= parameter(:σ_g, 2.5230, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_g: The standard deviation of the government spending process.",
                   tex_label="\\sigma_{g}")

    m <= parameter(:σ_b, 0.0292, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_b: The standard deviation of the intertemporal preference shifter process.",
                   tex_label="\\sigma_{b}")

    m <= parameter(:σ_μ, 0.4559, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
                   tex_label="\\sigma_{\\mu}")

    m <= parameter(:σ_ztil, 0.6742, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_ztil: The standard deviation of the process describing the stationary component of productivity.",
                   tex_label="\\sigma_{\\tilde{z}}")

    m <= parameter(:σ_λ_f, 0.1314, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_λ_f: The mean of the process that generates the price elasticity of the composite good. Specifically, the elasticity is (1+λ_{f,t})/(λ_{f_t}).",
                   tex_label="\\sigma_{\\lambda_f}")

    m <= parameter(:σ_λ_w, 0.3864, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{\\lambda_w}")

    m <= parameter(:σ_r_m, 0.2380, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_r_m: The standard deviation of the monetary policy shock.",
                   tex_label="\\sigma_{r^m}")

    m <= parameter(:σ_σ_ω, 0.0428, (0., 100.), (1e-5, 0.), ModelConstructors.Exponential(), RootInverseGamma(4, 0.05), fixed=false,
                   description="σ_σ_ω: The standard deviation of entrepreneurs' capital productivity follows an exogenous process with standard deviation σ_σ_ω.",
                   tex_label="\\sigma_{\\sigma_\\omega}")

    m <= parameter(:σ_μ_e, 0.0000, (0. ,100.), (1e-5, 0.), ModelConstructors.Exponential(), RootInverseGamma(4, 0.05), fixed=false,
                   description="σ_μ_e: Exogenous bankrupcy costs follow an exogenous process with standard deviation σ_μ_e.",
                   tex_label="\\sigma_{\\mu_e}")

    m <= parameter(:σ_γ, 0.0000, (0.,100.), (1e-5, 0.), ModelConstructors.Exponential(), RootInverseGamma(4, 0.01), fixed=false,
                   description="σ_γ: The fraction of entrepreneurs surviving period t follows an exogenous process with standard deviation σ_γ.",
                   tex_label="\\sigma_{\\gamma}")

    m <= parameter(:σ_π_star, 0.0269, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(6, 0.03), fixed=false,
                   description="σ_π_star: The standard deviation of the inflation target.",
                   tex_label="\\sigma_{\\pi_*}")

    m <= parameter(:σ_lr, 0.1766, (0.,10.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.75), fixed=false,
                   tex_label="\\sigma_{10y}")

    m <= parameter(:σ_z_p, 0.1662, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_z_p: The standard deviation of the shock to the permanent component of productivity.",
                   tex_label="\\sigma_{z^p}")

    m <= parameter(:σ_tfp, 0.9391, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{tfp}")

    m <= parameter(:σ_gdpdef, 0.1575, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{gdpdef}")

    m <= parameter(:σ_corepce, 0.0999, (0., 5.),(0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{pce}")

    m <= parameter(:σ_gdp, 0.1, (0., 5.),(0., 5.),ModelConstructors.Exponential(),RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{gdp}")

    m <= parameter(:σ_gdi, 0.1, (0., 5.),(0., 5.),ModelConstructors.Exponential(),RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{gdi}")

    # standard deviations of the anticipated policy shocks
    for i = 1:n_mon_anticipated_shocks(m)
        m <= parameter(Symbol("σ_r_m$(i)"), .2, (0., 100.), (1e-5, 0.), ModelConstructors.Exponential(),
                       RootInverseGamma(4, .2), fixed=false,
                       description="σ_r_m$(i): Standard deviation of the $i-period-ahead anticipated policy shock.",
                       tex_label=@sprintf("\\sigma_{ant%d}",i))
    end
    for i = n_mon_anticipated_shocks(m) + 1:n_mon_anticipated_shocks_padding(m)
        m <= parameter(Symbol("σ_r_m$(i)"), 0., (0., 100.), (1e-5, 0.), ModelConstructors.Exponential(),
                       RootInverseGamma(4, .2), fixed=false,
                       description="σ_r_m$(i): Standard deviation of the $i-period-ahead anticipated policy shock.",
                       tex_label=@sprintf("\\sigma_{ant%d}",i))
    end

    # Define regimes for standard shocks
    for i in 1:4
        ModelConstructors.set_regime_val!(m[:σ_g], i, m[:σ_g].value)
        ModelConstructors.set_regime_val!(m[:σ_b], i, m[:σ_b].value)
        ModelConstructors.set_regime_val!(m[:σ_μ], i, m[:σ_μ].value)
        ModelConstructors.set_regime_val!(m[:σ_ztil], i, m[:σ_ztil].value)
        ModelConstructors.set_regime_val!(m[:σ_λ_f], i, m[:σ_λ_f].value)
        ModelConstructors.set_regime_val!(m[:σ_λ_w], i, m[:σ_λ_w].value)
        ModelConstructors.set_regime_val!(m[:σ_r_m], i, m[:σ_r_m].value)
        ModelConstructors.set_regime_val!(m[:σ_σ_ω], i, m[:σ_σ_ω].value)
        ModelConstructors.set_regime_val!(m[:σ_μ_e], i, m[:σ_μ_e].value)
        ModelConstructors.set_regime_val!(m[:σ_γ], i, m[:σ_γ].value)
        ModelConstructors.set_regime_val!(m[:σ_π_star], i, m[:σ_π_star].value)
        ModelConstructors.set_regime_val!(m[:σ_lr], i, m[:σ_lr].value)
        ModelConstructors.set_regime_val!(m[:σ_z_p], i, m[:σ_z_p].value)
        ModelConstructors.set_regime_val!(m[:σ_tfp], i, m[:σ_tfp].value)
        err_adj = (i == 3) ? 10. : 1.
        ModelConstructors.set_regime_val!(m[:σ_gdpdef], i, err_adj * m[:σ_gdpdef].value)
        ModelConstructors.set_regime_val!(m[:σ_corepce], i, err_adj * m[:σ_corepce].value)
        ModelConstructors.set_regime_val!(m[:σ_gdp], i, m[:σ_gdp].value)
        ModelConstructors.set_regime_val!(m[:σ_gdi], i, m[:σ_gdi].value)

        for j = 1:DSGE.n_mon_anticipated_shocks(m)
            ModelConstructors.set_regime_val!(m[Symbol("σ_r_m$(j)")], i, m[Symbol("σ_r_m$(j)")].value)
        end
    end

    ModelConstructors.set_regime_val!(m[:σ_r_m], 2, 10. * m[:σ_r_m].value) # By default, we set these to 10 times their
    ModelConstructors.set_regime_val!(m[:σ_r_m], 3, 10. * m[:σ_r_m].value) # value in the pre-COVID and post-COVID regimes

    # Set up the COVID-19 shocks. We do not add anticipation by default
    ModelConstructors.set_regime_val!(m[:σ_φ], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_φ], 2, 400.) # initial value from literature, blown up by 100
    ModelConstructors.set_regime_val!(m[:σ_φ], 3, 400.)
    ModelConstructors.set_regime_val!(m[:σ_φ], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_φ_prop], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_φ_prop], 2, 2.)
    ModelConstructors.set_regime_val!(m[:σ_φ_prop], 3, 0.)
    ModelConstructors.set_regime_val!(m[:σ_φ_prop], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc], 2, 4.)
    ModelConstructors.set_regime_val!(m[:σ_biidc], 3, 4.)
    ModelConstructors.set_regime_val!(m[:σ_biidc], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc_prop], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc_prop], 2, 2.)
    ModelConstructors.set_regime_val!(m[:σ_biidc_prop], 3, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc_prop], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid], 2, 5.)
    ModelConstructors.set_regime_val!(m[:σ_ziid], 3, 5.)
    ModelConstructors.set_regime_val!(m[:σ_ziid], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid_prop], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid_prop], 2, 2.)
    ModelConstructors.set_regime_val!(m[:σ_ziid_prop], 3, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid_prop], 4, 0.)
end

"""
```
ss61!(m::Model1002)
```

models a "persistent demand shortfall" scenario of the COVID-19 pandemic by using the following features:

1. biidc, φ, and ziid shocks alive in 2020Q1-Q2
2. One-period ahead anticipated shock for all the iid shocks alive in 2020Q1 (expected to hit in 2020Q2). Only the
    one-period ahead anticipated shock for the biidc shock is alive in 2020Q2 (expected to hit in 2020Q3).
    The anticipated shocks in 2020Q1 are assumed to be twice the size of
    the unanticipated shocks realized in 2020Q1. The anticipated biidc shock in 2020Q2 is assumed to be the same size of
    the unanticipated biidc shock realized in 2020Q2.
3. The standard deviation of monetary policy shocks are set to 10 times their value in the time periods before 2020Q1.

Confer with the official documentation about the COVID-19 shocks on the GitHub page.
"""
function ss61!(m::Model1002)
    # Set up regime-switching
    m <= Setting(:regime_switching, true)
    m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2020, 3, 31),
                                                3 => Date(2020, 6, 30), 4 => Date(2020, 9, 30)))
    setup_regime_switching_inds!(m)

    # Allow the lower bound of non-COVID-19 parameters to equal zero
    m <= parameter(:σ_g, 2.5230, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_g: The standard deviation of the government spending process.",
                   tex_label="\\sigma_{g}")

    m <= parameter(:σ_b, 0.0292, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_b: The standard deviation of the intertemporal preference shifter process.",
                   tex_label="\\sigma_{b}")

    m <= parameter(:σ_μ, 0.4559, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
                   tex_label="\\sigma_{\\mu}")

    m <= parameter(:σ_ztil, 0.6742, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_ztil: The standard deviation of the process describing the stationary component of productivity.",
                   tex_label="\\sigma_{\\tilde{z}}")

    m <= parameter(:σ_λ_f, 0.1314, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_λ_f: The mean of the process that generates the price elasticity of the composite good. Specifically, the elasticity is (1+λ_{f,t})/(λ_{f_t}).",
                   tex_label="\\sigma_{\\lambda_f}")

    m <= parameter(:σ_λ_w, 0.3864, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{\\lambda_w}")

    m <= parameter(:σ_r_m, 0.2380, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_r_m: The standard deviation of the monetary policy shock.",
                   tex_label="\\sigma_{r^m}")

    m <= parameter(:σ_σ_ω, 0.0428, (0., 100.), (1e-5, 0.), ModelConstructors.Exponential(), RootInverseGamma(4, 0.05), fixed=false,
                   description="σ_σ_ω: The standard deviation of entrepreneurs' capital productivity follows an exogenous process with standard deviation σ_σ_ω.",
                   tex_label="\\sigma_{\\sigma_\\omega}")

    m <= parameter(:σ_μ_e, 0.0000, (0. ,100.), (1e-5, 0.), ModelConstructors.Exponential(), RootInverseGamma(4, 0.05), fixed=false,
                   description="σ_μ_e: Exogenous bankrupcy costs follow an exogenous process with standard deviation σ_μ_e.",
                   tex_label="\\sigma_{\\mu_e}")

    m <= parameter(:σ_γ, 0.0000, (0.,100.), (1e-5, 0.), ModelConstructors.Exponential(), RootInverseGamma(4, 0.01), fixed=false,
                   description="σ_γ: The fraction of entrepreneurs surviving period t follows an exogenous process with standard deviation σ_γ.",
                   tex_label="\\sigma_{\\gamma}")

    m <= parameter(:σ_π_star, 0.0269, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(6, 0.03), fixed=false,
                   description="σ_π_star: The standard deviation of the inflation target.",
                   tex_label="\\sigma_{\\pi_*}")

    m <= parameter(:σ_lr, 0.1766, (0.,10.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.75), fixed=false,
                   tex_label="\\sigma_{10y}")

    m <= parameter(:σ_z_p, 0.1662, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_z_p: The standard deviation of the shock to the permanent component of productivity.",
                   tex_label="\\sigma_{z^p}")

    m <= parameter(:σ_tfp, 0.9391, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{tfp}")

    m <= parameter(:σ_gdpdef, 0.1575, (0., 5.), (0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{gdpdef}")

    m <= parameter(:σ_corepce, 0.0999, (0., 5.),(0., 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{pce}")

    m <= parameter(:σ_gdp, 0.1, (0., 5.),(0., 5.),ModelConstructors.Exponential(),RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{gdp}")

    m <= parameter(:σ_gdi, 0.1, (0., 5.),(0., 5.),ModelConstructors.Exponential(),RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{gdi}")

    # standard deviations of the anticipated policy shocks
    for i = 1:n_mon_anticipated_shocks(m)
        m <= parameter(Symbol("σ_r_m$(i)"), .2, (0., 100.), (1e-5, 0.), ModelConstructors.Exponential(),
                       RootInverseGamma(4, .2), fixed=false,
                       description="σ_r_m$(i): Standard deviation of the $i-period-ahead anticipated policy shock.",
                       tex_label=@sprintf("\\sigma_{ant%d}",i))
    end
    for i = n_mon_anticipated_shocks(m) + 1:n_mon_anticipated_shocks_padding(m)
        m <= parameter(Symbol("σ_r_m$(i)"), 0., (0., 100.), (1e-5, 0.), ModelConstructors.Exponential(),
                       RootInverseGamma(4, .2), fixed=false,
                       description="σ_r_m$(i): Standard deviation of the $i-period-ahead anticipated policy shock.",
                       tex_label=@sprintf("\\sigma_{ant%d}",i))
    end

    # Define regimes for standard shocks
    for i in 1:4
        ModelConstructors.set_regime_val!(m[:σ_g], i, m[:σ_g].value)
        ModelConstructors.set_regime_val!(m[:σ_b], i, m[:σ_b].value)
        ModelConstructors.set_regime_val!(m[:σ_μ], i, m[:σ_μ].value)
        ModelConstructors.set_regime_val!(m[:σ_ztil], i, m[:σ_ztil].value)
        ModelConstructors.set_regime_val!(m[:σ_λ_f], i, m[:σ_λ_f].value)
        ModelConstructors.set_regime_val!(m[:σ_λ_w], i, m[:σ_λ_w].value)
        ModelConstructors.set_regime_val!(m[:σ_r_m], i, m[:σ_r_m].value)
        ModelConstructors.set_regime_val!(m[:σ_σ_ω], i, m[:σ_σ_ω].value)
        ModelConstructors.set_regime_val!(m[:σ_μ_e], i, m[:σ_μ_e].value)
        ModelConstructors.set_regime_val!(m[:σ_γ], i, m[:σ_γ].value)
        ModelConstructors.set_regime_val!(m[:σ_π_star], i, m[:σ_π_star].value)
        ModelConstructors.set_regime_val!(m[:σ_lr], i, m[:σ_lr].value)
        ModelConstructors.set_regime_val!(m[:σ_z_p], i, m[:σ_z_p].value)
        ModelConstructors.set_regime_val!(m[:σ_tfp], i, m[:σ_tfp].value)
        err_adj = (i == 3) ? 10. : 1.
        ModelConstructors.set_regime_val!(m[:σ_gdpdef], i, err_adj * m[:σ_gdpdef].value)
        ModelConstructors.set_regime_val!(m[:σ_corepce], i, err_adj * m[:σ_corepce].value)
        ModelConstructors.set_regime_val!(m[:σ_gdp], i, m[:σ_gdp].value)
        ModelConstructors.set_regime_val!(m[:σ_gdi], i, m[:σ_gdi].value)

        for j = 1:DSGE.n_mon_anticipated_shocks(m)
            ModelConstructors.set_regime_val!(m[Symbol("σ_r_m$(j)")], i, m[Symbol("σ_r_m$(j)")].value)
        end
    end

    ModelConstructors.set_regime_val!(m[:σ_r_m], 2, 10. * m[:σ_r_m].value) # By default, we set these to 10 times their
    ModelConstructors.set_regime_val!(m[:σ_r_m], 3, 10. * m[:σ_r_m].value) # value in the pre-COVID and post-COVID regimes

    # Set up the COVID-19 shocks. We do not add anticipation by default
    ModelConstructors.set_regime_val!(m[:σ_φ], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_φ], 2, 400.) # initial value from literature, blown up by 100
    ModelConstructors.set_regime_val!(m[:σ_φ], 3, 400.)
    ModelConstructors.set_regime_val!(m[:σ_φ], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_φ_prop], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_φ_prop], 2, 2.)
    ModelConstructors.set_regime_val!(m[:σ_φ_prop], 3, 0.)
    ModelConstructors.set_regime_val!(m[:σ_φ_prop], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc], 2, 4.)
    ModelConstructors.set_regime_val!(m[:σ_biidc], 3, 4.)
    ModelConstructors.set_regime_val!(m[:σ_biidc], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc_prop], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc_prop], 2, 2.)
    ModelConstructors.set_regime_val!(m[:σ_biidc_prop], 3, 1.)
    ModelConstructors.set_regime_val!(m[:σ_biidc_prop], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid], 2, 5.)
    ModelConstructors.set_regime_val!(m[:σ_ziid], 3, 5.)
    ModelConstructors.set_regime_val!(m[:σ_ziid], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid_prop], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid_prop], 2, 2.)
    ModelConstructors.set_regime_val!(m[:σ_ziid_prop], 3, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid_prop], 4, 0.)
end

ext: jl


text:
"""
   KnownDerivativesMonotonicInterpolation(m)

Cubic interpolations with given derivatives at nodes.
"""
struct KnownDerivativesMonotonicInterpolation{TEl} <:
       Interpolations.MonotonicInterpolationType
    m::Vector{TEl}
end

function Interpolations.calcTangents(
    ::Type{TCoeffs},
    x::AbstractVector{<:Number},
    y::AbstractVector{TEl},
    method::KnownDerivativesMonotonicInterpolation,
) where {TCoeffs,TEl}
    n = length(x)
    Δ = Vector{TCoeffs}(undef, n - 1)
    m = Vector{TCoeffs}(undef, n)
    for k in 1:(n - 1)
        Δ[k] = (y[k + 1] - y[k]) / (x[k + 1] - x[k])
        m[k] = method.m[k]
    end
    m[n] = method.m[n]
    return (m, Δ)
end

function ManifoldsBase.number_eltype(itp::Interpolations.AbstractInterpolation)
    return eltype(itp)
end


text:
function dot_simd(a::AbstractVector, b::AbstractVector)
    s = zero(eltype(a))
    @fastmath @inbounds @simd for i ∈ eachindex(a)
        s += a[i]' * b[i]
    end
    s
end
function cdot_mat(ca::AbstractVector{Complex{T}}, cb::AbstractVector{Complex{T}}) where {T}
    a = reinterpret(reshape, T, ca)
    b = reinterpret(reshape, T, cb)
    re = zero(T); im = zero(T)
    @avx for i ∈ axes(a,2)
        re += a[1,i] * b[1,i] + a[2,i] * b[2,i]
        im += a[1,i] * b[2,i] - a[2,i] * b[1,i]
    end
    Complex(re, im)
end
function cdot_affine(ca::AbstractVector{Complex{T}}, cb::AbstractVector{Complex{T}}) where {T}
    a = reinterpret(T, ca);
    b = reinterpret(T, cb);
    re = zero(T); im = zero(T)
    # with a multiplier, we go from `i = 1 -> 2i = 2` to `i = 0 -> 2i = 0
    # 2(i+1-1) = 2i + 2 - 2, so....
    @avx for i ∈ 1:length(a)>>>1
        re += a[2i-1] * b[2i-1] + a[2i] * b[2i  ]
        im += a[2i-1] * b[2i  ] - a[2i] * b[2i-1]
    end
    Complex(re, im)
end
function cdot_stride(ca::AbstractVector{Complex{T}}, cb::AbstractVector{Complex{T}}) where {T}
    a = reinterpret(T, ca);
    b = reinterpret(T, cb);
    re = zero(T); im = zero(T)
    @avx for i ∈ 1:2:length(a)
        re += a[i] * b[i  ] + a[i+1] * b[i+1]
        im += a[i] * b[i+1] - a[i+1] * b[i  ]
    end
    Complex(re, im)
end
function qdot_simd(x::AbstractVector{NTuple{4,T}}, y::AbstractVector{NTuple{4,T}}) where {T}
    a = zero(T)
    b = zero(T)
    c = zero(T)
    d = zero(T)
    @fastmath @inbounds @simd for i ∈ eachindex(x)
        a₁, b₁, c₁, d₁ = x[i]
        a₂, b₂, c₂, d₂ = y[i]
        a += a₁*a₂ + b₁*b₂ + c₁*c₂ + d₁*d₂
        b += a₁*b₂ - b₁*a₂ - c₁*d₂ + d₁*c₂
        c += a₁*c₂ + b₁*d₂ - c₁*a₂ - d₁*b₂
        d += a₁*d₂ - b₁*c₂ + c₁*b₂ - d₁*a₂
    end
    (a,b,c,d)
end
function qdot_mat(x::AbstractMatrix, y::AbstractMatrix)
    a = zero(eltype(x))
    b = zero(eltype(x))
    c = zero(eltype(x))
    d = zero(eltype(x))
    @avx for i ∈ axes(x,2)
        a₁ = x[1,i]
        b₁ = x[2,i]
        c₁ = x[3,i]
        d₁ = x[4,i]
        a₂ = y[1,i]
        b₂ = y[2,i]
        c₂ = y[3,i]
        d₂ = y[4,i]
        a += a₁*a₂ + b₁*b₂ + c₁*c₂ + d₁*d₂
        b += a₁*b₂ - b₁*a₂ - c₁*d₂ + d₁*c₂
        c += a₁*c₂ + b₁*d₂ - c₁*a₂ - d₁*b₂
        d += a₁*d₂ - b₁*c₂ + c₁*b₂ - d₁*a₂
    end
    (a,b,c,d)
end
function qdot_affine(x::AbstractVector, y::AbstractVector)
    a = zero(eltype(x))
    b = zero(eltype(x))
    c = zero(eltype(x))
    d = zero(eltype(x))
    @avx for i ∈ 1:length(x)>>2
        a₁ = x[4i-3]
        b₁ = x[4i-2]
        c₁ = x[4i-1]
        d₁ = x[4i]
        a₂ = y[4i-3]
        b₂ = y[4i-2]
        c₂ = y[4i-1]
        d₂ = y[4i]
        a += a₁*a₂ + b₁*b₂ + c₁*c₂ + d₁*d₂
        b += a₁*b₂ - b₁*a₂ - c₁*d₂ + d₁*c₂
        c += a₁*c₂ + b₁*d₂ - c₁*a₂ - d₁*b₂
        d += a₁*d₂ - b₁*c₂ + c₁*b₂ - d₁*a₂
    end
    (a,b,c,d)
end
function qdot_stride(x::AbstractVector, y::AbstractVector)
    a = zero(eltype(x))
    b = zero(eltype(x))
    c = zero(eltype(x))
    d = zero(eltype(x))
    @avx for i ∈ 1:4:length(x)
        a₁ = x[i]
        b₁ = x[i+1]
        c₁ = x[i+2]
        d₁ = x[i+3]
        a₂ = y[i]
        b₂ = y[i+1]
        c₂ = y[i+2]
        d₂ = y[i+3]
        a += a₁*a₂ + b₁*b₂ + c₁*c₂ + d₁*d₂
        b += a₁*b₂ - b₁*a₂ - c₁*d₂ + d₁*c₂
        c += a₁*c₂ + b₁*d₂ - c₁*a₂ - d₁*b₂
        d += a₁*d₂ - b₁*c₂ + c₁*b₂ - d₁*a₂
    end
    (a,b,c,d)
end
function cmatmul_array!(C::AbstractArray{T,3}, A::AbstractArray{T,3}, B::AbstractArray{T,3}) where {T}
    @avx for n ∈ indices((C,B),3), m ∈ indices((C,A),2)
        Cre = zero(T)
        Cim = zero(T)
        for k ∈ indices((A,B),(3,2))
            Cre += A[1,m,k] * B[1,k,n] - A[2,m,k] * B[2,k,n]
            Cim += A[1,m,k] * B[2,k,n] + A[2,m,k] * B[1,k,n]
        end
        C[1,m,n] = Cre
        C[2,m,n] = Cim
    end
end

function issue209(M, G, J, H, A, B, ϕ)
    tmp = similar(ϕ, G-1, (2*J+1)*(H + 1))
    Bf = reinterpret(reshape, Float64, B)
    ϕf = reinterpret(reshape, Float64, ϕ)
    tmpf = reinterpret(reshape, Float64, tmp)
    jmax = 2*J + 1
    for mm = 1:M
        m_idx = M + 2 - mm
        @avx for hh = 1:H+1
            h_idx = (hh - 1)*jmax
            for jj = 1:jmax, gg = 1:G-1
                tmpf[1, gg, jj + h_idx] = ϕf[1, jj, gg+1, hh, m_idx] +
                                          Bf[1, jj, gg,   hh, m_idx]
                tmpf[2, gg, jj + h_idx] = ϕf[2, jj, gg+1, hh, m_idx] +
                                          Bf[2, jj, gg,   hh, m_idx]
            end
        end
    end
    tmp
end
function issue209_noavx(M, G, J, H, A, B, ϕ)
    tmp = similar(ϕ, G-1, (2*J+1)*(H + 1))
    Bf = reinterpret(reshape, Float64, B)
    ϕf = reinterpret(reshape, Float64, ϕ)
    tmpf = reinterpret(reshape, Float64, tmp)
    jmax = 2*J + 1
    for mm = 1:M
        m_idx = M + 2 - mm
        for hh = 1:H+1
            h_idx = (hh - 1)*jmax
            for jj = 1:jmax, gg = 1:G-1
                tmpf[1, gg, jj + h_idx] = ϕf[1, jj, gg+1, hh, m_idx] +
                                          Bf[1, jj, gg,   hh, m_idx]
                tmpf[2, gg, jj + h_idx] = ϕf[2, jj, gg+1, hh, m_idx] +
                                          Bf[2, jj, gg,   hh, m_idx]
            end
        end
    end
    tmp
end

@testset "shuffles load/stores" begin
	@show @__LINE__
    for i ∈ 1:128
        ac = rand(Complex{Float64}, i);
        bc = rand(Complex{Float64}, i);
        dsimd = dot_simd(ac, bc)
        if VERSION ≥ v"1.6.0-rc1"
            @test dsimd ≈ cdot_mat(ac, bc)
        end
        @test dsimd ≈ cdot_affine(ac, bc) ≈ cdot_stride(ac, bc)

        xq = [ntuple(_ -> rand(), Val(4)) for _ ∈ 1:i];
        yq = [ntuple(_ -> rand(), Val(4)) for _ ∈ 1:i];
        xqv = reinterpret(Float64, xq);
        yqv = reinterpret(Float64, yq);
        qsimd = Base.vect(qdot_simd(xq, yq)...);
        if VERSION ≥ v"1.6.0-rc1"
            xqm = reinterpret(reshape, Float64, xq);
            yqm = reinterpret(reshape, Float64, yq);
            @test qsimd ≈ Base.vect(qdot_mat(xqm, yqm)...)
        end
        @test qsimd ≈ Base.vect(qdot_affine(xqv, yqv)...) ≈ Base.vect(qdot_stride(xqv, yqv)...)

        if VERSION ≥ v"1.6.0-rc1"
            Ac = rand(Complex{Float64}, i, i);
            Bc = rand(Complex{Float64}, i, i);
            Cc1 = Ac*Bc;
            Cc2 = similar(Cc1);
            # Cc3 = similar(Cc1)
            Cca = reinterpret(reshape, Float64, Cc2);
            Aca = reinterpret(reshape, Float64, Ac);
            Bca = reinterpret(reshape, Float64, Bc);
            cmatmul_array!(Cca, Aca, Bca)

            @test Cc1 ≈ Cc2# ≈ Cc3
        end
    end
    if VERSION ≥ v"1.6.0-rc1"
        M = 100
        G = 50
        J = 50
        H = 300

        A = Matrix(Tridiagonal(rand(G-1,G-1)));
        B = rand(Complex{Float64}, 2*J+1, G-1, H+1, M+1);
        ϕ = rand(Complex{Float64}, 2*J+1, G+1, H+1, M+1);
        @test issue209(M, G, J, H, A, B, ϕ) ≈ issue209_noavx(M, G, J, H, A, B, ϕ)
    end
end


text:
#!/usr/bin/env julia

using Pkg
Pkg.activate(".")
#Pkg.instantiate()
Pkg.status()

time_start = time()
println("start pre-compilation")

using Distributed
using SparseArrays

using JSON

using JuMP

using Ipopt

using InfrastructureModels
using PowerModels
using Memento

println("pre-compilation finished ($(time() - time_start))")


text:
using FFTW
using LoopVectorization
using MCMCDiagnosticTools

using Tullio

export relative_eff, psis_ess, sup_ess

"""
    relative_eff(
        sample::AbstractArray{Real, 3};
        method=MCMCDiagnosticTools.FFTESSMethod()
    )

Calculate the relative efficiency of an MCMC chain, i.e. the effective sample size divided
by the nominal sample size.
"""
function relative_eff(
    sample::AbstractArray{T, 3}; method=MCMCDiagnosticTools.FFTESSMethod()
) where {T <: Union{Real, Missing}}
    dims = size(sample)
    post_sample_size = dims[2] * dims[3]
    ess_sample = inv.(permutedims(sample, [2, 1, 3]))
    ess, = MCMCDiagnosticTools.ess_rhat(ess_sample; method=method, maxlag=dims[2])
    r_eff = ess / post_sample_size
    return r_eff
end

"""
    function psis_ess(
        weights::AbstractVector{T},
        r_eff::AbstractVector{T}
    ) -> AbstractVector{T}

Calculate the (approximate) effective sample size of a PSIS sample, using the correction in
Vehtari et al. 2019. This uses the variance-based definition of ESS, and measures the L2
distance of the proposal and target distributions.

# Arguments

  - `weights`: A set of importance sampling weights derived from PSIS.
  - `r_eff`: The relative efficiency of the MCMC chains from which PSIS samples were derived.

See `?relative_eff` to calculate `r_eff`.
"""
function psis_ess(
    weights::AbstractVector{T}, r_eff::AbstractVector{T}
) where {T <: Union{Real, Missing}}
    @tullio sum_of_squares := weights[x]^2
    return @turbo r_eff ./ sum_of_squares
end

function psis_ess(
    weights::AbstractMatrix{T}, r_eff::AbstractVector{T}
) where {T <: Union{Real, Missing}}
    @tullio sum_of_squares[x] := weights[x, y]^2
    return @turbo r_eff ./ sum_of_squares
end

function psis_ess(weights::AbstractMatrix{<:Union{Real, Missing}})
    @warn "PSIS ESS not adjusted based on MCMC ESS. MCSE and ESS estimates " *
          "will be overoptimistic if samples are autocorrelated."
    return psis_ess(weights, ones(size(weights)))
end

"""
    function sup_ess(
        weights::AbstractVector{T},
        r_eff::AbstractVector{T}
    ) -> AbstractVector

Calculate the supremum-based effective sample size of a PSIS sample, i.e. the inverse of the
maximum weight. This measure is more trustworthy than the `ess` from `psis_ess`. It uses the
L-∞ norm.

# Arguments
  - `weights`: A set of importance sampling weights derived from PSIS.
  - `r_eff`: The relative efficiency of the MCMC chains from which PSIS samples were derived.
"""
function sup_ess(
    weights::AbstractMatrix{T}, r_eff::V
 ) where {T<:Union{Real, Missing}, V<:AbstractVector{T}}
    return @turbo inv.(dropdims(maximum(weights; dims=2); dims=2)) .* r_eff
end


text:
#export AbstractSymmetryOperation
using GroupTools

export AbstractSpaceSymmetryOperation
export domaintype

# abstract type AbstractSymmetryOperation end

"""
    AbstractSpaceSymmetryOperation{S<:Real}

Abstract space symmetry operation, i.e. translation, point, and space operation.
See also [`TranslationOperation`](@ref), [`PointOperation`](@ref), [`SpaceOperation`](@ref).
"""
abstract type AbstractSpaceSymmetryOperation{S<:Real} <: AbstractSymmetryOperation end

"""
    domaintype(arg::AbstractSpaceSymmetryOperation{S}) where {S<:Real}

Domain type of `arg`, i.e. the type of the coordinates.
"""
domaintype(::AbstractSpaceSymmetryOperation{S}) where {S<:Real} = S
domaintype(::Type{<:AbstractSpaceSymmetryOperation{S}}) where {S<:Real} = S


text:
using Dagger
import Dagger: reduceblock
using DistributedArrays

using BenchmarkTools

AmulB_compatible(x,y) =
    try; size(x), size(y), x * y; true; catch err false end

AmulB_compatible(x::Computed,y::Computed) = AmulB_compatible(domain(x.result), domain(y.result))

function AcmulB_compatible(x,y)
    AmulB_compatible(domain(x.result)', domain(y.result))
end

function matrix_suite(ctx, X, f=x->x)
    suite = BenchmarkGroup(["matrix"])

    suite["alloc"] = @benchmarkable (compute($(f(X))); gc())

    Y = compute(ctx, f(X))
    suite["X+1.0"]  = @benchmarkable (compute($ctx, $(f(Y+1))); gc())
    suite["X+X"]    = @benchmarkable (compute($ctx, $(f(Y+Y))); gc())
    suite["X'"]     = @benchmarkable (compute($ctx, $(f(Y'))); gc())
    suite["sin(X)"] = @benchmarkable (compute($ctx, $(f(sin(Y)))); gc())

    if domain(Y.result) == domain(Y.result)'
        @show domain(Y.result)
        suite["X+X'"]   = @benchmarkable (compute($ctx, $(f(Y+Y'))); gc())
    end
    if AmulB_compatible(Y, Y)
        suite["X*X"]    = @benchmarkable (compute($ctx, $(f(Y*Y))); gc())
    end
    if AcmulB_compatible(Y, Y)
        suite["X'*X"]   = @benchmarkable (compute($ctx, $(f(Y'*Y))); gc())
    end
    suite
end

function matrix_variants(NW, f=x->x)
    w = workers()[1:NW]
    suite = BenchmarkGroup()

    matrix(ch, sz) = zeros(Chunks(ch), sz)

    suite["trivial_one_each"] = matrix_suite(Context(w), matrix(1, NW), f)
    suite["trivial_two_each"] = matrix_suite(Context(w), matrix(1, 2NW), f)
    suite["trivial_1k_each"] = matrix_suite(Context(w), matrix(1, 1000NW), f)

 #  suite["128M_one_each"] =
 #      matrix_suite(Context(w), matrix((4000,4000), (4000*NW, 4000)), f)

 #  suite["128M_two_each"] =
 #      matrix_suite(Context(w), matrix((4000,4000), (NW*4000, 2*4000)), f)

 #  suite["128M_uneven"] =
 #      matrix_suite(Context(w), matrix((4000,4000), (5000*NW, 5000)), f)
    suite
end

function matrix_suite_nworkers(NW)
    suite = BenchmarkGroup()
    suite["in-memory"] = matrix_variants(2, x->x)
    suite["on-disk"]   = matrix_variants(2, x -> Dagger.save(x, "tmp$(randstring(5))"))
    suite
end

using JLD
function main()
    if length(ARGS) == 0
        error("Usage: julia -p NWorkers benchmark.jl [tune-file] [output-file]")
    end

    nw = nworkers()
    suite = matrix_suite_nworkers(nw)
    if length(ARGS) > 0
        tune_file = ARGS[1]
        if isfile(tune_file)
            println("Using benchmark tuning data in $tune_file")
            loadparams!(suite, JLD.load(tune_file, "suite"), :evals, :samples)
        else
            println("Creating benchmark tuning file $tune_file")
            tune!(suite)
            JLD.save(tune_file, "suite", params(suite))
        end
    else
        println("Tuning benchmarks. Tip: provide a tune_file argument to save this step.")
        tune!(suite)
    end

    println("running benchmarks")
    res = run(suite)
    @show(res)
    output = get(ARGS, 2, "result-$(nw)worker-$(Dates.now()).jld")
    println("saving results in $output")
    JLD.save(output, "result", res)
    println("Done.")

    if length(ARGS) > 2
        compare_file = ARGS[3]
        println("Comparing with results in $compare_file")
        reference_results = JLD.load(compare_file, "result")
        @show judge(reference_results, res)
    end

end

main()


text:
#######################
# This accumulates a stochastic estimate
# of the score using the algorithm of
# Nemeth et al. (2017).
export score, score_filter

struct ScoreParticleContainer <: ParticleContainer
    X # Particles
    w # weights
    α
    ℓ # Estimate of the log likelihood
    S # Score
    dg # Preallocated GradientResult
end

Base.size(ps::ScoreParticleContainer) = (size(ps.X,2),)
Base.length(ps::ScoreParticleContainer) = prod(size(ps))
Base.getindex(ps::ScoreParticleContainer,i::Int) = ps.X[:,i]
Base.setindex!(ps::ScoreParticleContainer,v,i::Int) = ps.X[:,i] = v
Base.IndexStyle(ps::ScoreParticleContainer) = IndexLinear()

StatsBase.score(ps::ScoreParticleContainer) = ps.S
StatsBase.loglikelihood(ps::ScoreParticleContainer) = ps.ℓ[end]

resample!(ps::ScoreParticleContainer,m,t,θ,y,u,sample_function=resample_stratified) = resample!(ScoreParticleContainer,ps.X,ps.w,ps.α,ps.w[:,t],m,t,θ,y,u,sample_function)

function resample!(::Type{ScoreParticleContainer},X,w,α,ξs,m,t,θ,y,u,sample_function)
    for i in eachindex(ξs)
        ξs[i] += weight_function(m,X[:,i,t],y,u,t,θ)
    end
    Σ = sum(exp,ξs)
    idx = sample_function(exp.(ξs)./Σ,length(ξs))

    X[:,:,t] .= X[:,idx,t]
    w[:,t] .= w[idx,t]
    α[:,:,t] .= α[:,idx,t]
    ξs[idx],Σ
end

propagate!(ps::ScoreParticleContainer,ξ,Σ,m,t,θ,y,u,λ) = propagate!(ScoreParticleContainer,ps.X,ps.w,ps.α,ps.ℓ,ps.S,ps.dg,ξ,Σ,m,t,θ,y,u,λ)

function propagate!(::Type{ScoreParticleContainer},X,w,α,ℓ,S,dg,ξ,Σ,m,t,θ,y,u,λ)
    for i in eachindex(ξ)
        #dq = proposal_distribution(m,X[:,i,t],y,u,t,θ)
        #rand!(dq,view(X,:,i,t+1))
        proposal_rand!(view(X,:,i,t+1),m,X[:,i,t],y,u,t,θ)

        # We will need to define this more clearly
        dg = ssm_gradient!(dg,m,X[:,i,t+1],X[:,i,t],y,u,t,θ)

        w[i,t+1] = w[i,t]
        # _ssm_gradient!(...) should also calculate the
        # value of the two logpdfs.
        w[i,t+1] += DiffResults.value(dg)
        #w[i,t+1] += transition_logpdf(m,X[:,i,t+1],X[:,i,t],u,t,θ)
        #w[i,t+1] += observation_logpdf(m,y,X[:,i,t+1],u,t,θ)
        #w[i] *= pdf(f(X[:,i],t,θ),Xn)
        #w[i] *= pdf(g(Xn,t,θ),y)
        w[i,t+1] -= ξ[i]
        w[i,t+1] -= proposal_logpdf(m,X[:,i,t+1],X[:,i,t],y,u,t,θ)

        # This part is unique to the score filter
        α[:,i,t+1] .= α[:,i,t]
        α[:,i,t+1] .*= λ
        α[:,i,t+1] .+= (1-λ)*S[:] .+ DiffResults.gradient(dg)
    end
    ℓn = log(sum(exp,w[:,t+1]))
    ℓ[t+1] = ℓ[t] + ℓn - log(length(ξ)) + log(Σ)
    w[:,t+1] .-= ℓn
    S[:] .= (α[:,:,t+1]*exp.(w[:,t+1]))
end

function step(ps::ScoreParticleContainer,m::StateSpaceModel,t,θ,y,u;sample_function=resample_stratified,λ=0.95)

    # Resample

    #rs = Resampler(ps,m,y,t,θ,sample_function)

    ξs,Σ = resample!(ps,m,t,θ,y,u,sample_function)

    # Propagate

    #ps2 = ScoreParticleContainer([propose(t,x,w,ξ,α,ps.S,y,θ,m.f,m.g,m.q,m.dp,ps.dg,λ) for (x,w,ξ,α) in rs],ps.ℓ,ps.dg,sum(rs.ξ),λ)

    propagate!(ps,ξs,Σ,m,t,θ,y,u,λ)

    #normalize_weights!(ps2)

    ps
end

function score_filter(m::StateSpaceModel,θ,Y,N,u=fill([0.0],length(Y));sample_function=resample_stratified,verbose=false,λ=0.95)

    T = length(Y)

    if verbose
        p = Progress(T,1)
    end

    ps = initialize_particles(ScoreParticleContainer,m,θ,Y,N)

    for t in 1:T
        ps = step(ps,m,t,θ,Y[t],u[t],sample_function=sample_function,λ=λ)
        if verbose
            next!(p)
        end
    end

    ps
end

function initialize_particles(::Type{ScoreParticleContainer},m::StateSpaceModel,θ,Y,N)

    # We should probably think of μ0 as a distribution for p(x_0)
    # not one for p(x_1)
    # Initialization gives X0, w = 1/N
    #x0 = rand(μ0(θ),N)
    T = length(Y)

    x0 = initial_rand(m,θ,N)
    X0 = zeros(size(x0,1),N,T+1)
    X0[:,:,1] = x0[:,:]

    w = log.(ones(N,T+1)/N)
    ℓ = zeros(T+1)
    #w = [pdf(g(X0[:,end,i],1,θ),y)/N for i in 1:N]
    #ℓ = sum(w)
    #w ./= ℓ

    dg = DiffResults.GradientResult(θ)

    α = Array{Float64}(undef,length(θ),N,T+1)

    for i in 1:N
        # TODO: Implement this
        dg = initial_gradient!(dg,m,X0[:,i,1],θ)

        α[:,i,1] = DiffResults.gradient(dg)
    end

    S = α[:,:,1] * exp.(w[:,1])

    ScoreParticleContainer(X0,w,α,ℓ,S,dg)
end

#=
# Does it make sense to FB smooth the score filter?
# We can just make it identical to the particle smoother
# or we can figure out how to smooth the gradients as well

function particle_smoother(m::StateSpaceModel,θ,Y,pf::ScoreParticleContainer,B,u=fill([0.0],length(Y));sample_function=resample_stratified)
    D = size(pf.X,1)
    N = size(pf.X,2)
    T = size(pf.X,3)-1

    x = zeros(D,T+1,B)
    for i in 1:B
        x[:,end,i] = pf.X[:,sample_function(exp.(pf.w[:,end])./sum(exp,pf.w[:,end]),1),end]
    end

    for t in T-1:-1:0
        wsmooth = pf.w[:,t+1]
        for j in 1:B
            for i in 1:N
            # This is quite slow
            # We still need u!
                wsmooth[i] += logpdf(transition_distribution(m,pf.X[:,i,t+1],u[t+1],t,θ),x[:,t+2,j])
            end

            x[:,t+1,j] = pf.X[:,sample_function(exp.(wsmooth)./sum(exp,wsmooth),1),t+1]
        end
    end
    x
end
=#


text:
Base.run(`clear`)
println("Starting")
println("----------------------------")
using Leapfrog
println("JIT compilation")
@elapsed 1
mp_burn, cost_burn = init_all(1)
solve(cost_burn, mp_burn)
mp, cost = init_all(2);
@elapsed stags, ess, number_of_equilibria = solve(cost, mp)
println("----------------------------")
println("Checking # of equilibria")
println(number_of_equilibria)
println("----------------------------")


text:
################################
# Complex Operator Definitions #
################################
@present ExtendedOperators <: ExtendedFreeExtCalc2D begin
  X::Space
  neg::Hom(DualForm1(X), DualForm1(X)) # negative
  neg₁::Hom(Form1(X), Form1(X)) # negates 1-forms
  dneg₁::Hom(DualForm1(X), DualForm1(X)) # negates 1-forms
  L₀::Hom(Form1(X)⊗DualForm2(X), DualForm2(X))
  L₁::Hom(Form1(X)⊗DualForm1(X), DualForm1(X))
  i₀::Hom(Form1(X)⊗DualForm2(X), DualForm1(X))
  i₁::Hom(Form1(X)⊗DualForm1(X), DualForm0(X))

  L₁′::Hom(Form1(X)⊗Form1(X), Form1(X))
  i₀′::Hom(Form1(X)⊗Form2(X), Form1(X))
  i₁′::Hom(Form1(X)⊗Form1(X), Form0(X))
  ∧₁₁′::Hom(Form1(X)⊗DualForm1(X), DualForm2(X))
  ∧₁₀′::Hom(Form1(X)⊗DualForm0(X), Form1(X))
end

rules = gen_dec_rules()

# Lie derivative between two primal 1-forms
Lie1Imp′ = @decapode ExtendedOperators begin
  (F1, F1′, Ḟ1)::Form1{X}
  Ḟ1 == i₀′(F1, d₁{X}(F1′)) + d₀{X}(i₁′(F1, F1′))
end
lie1_imp′ = diag2dwd(Lie1Imp′, in_vars = [:F1, :F1′], out_vars = [:Ḟ1])
rules[:L₁′] = lie1_imp′

# Internal product between a primal 1-form and a primal 2-form
I0Imp′ = @decapode ExtendedOperators begin
  (F1, Ḟ1)::Form1{X}
  F2::Form2{X}
  Ḟ1 == neg₁(∧₁₀′(F1, ⋆₂{X}(F2)))
end
i0_imp′ = diag2dwd(I0Imp′, in_vars = [:F1, :F2], out_vars = [:Ḟ1])
rules[:i₀′] = i0_imp′

# Internal product between two primal 1-forms
I1Imp′ = @decapode ExtendedOperators begin
  (F1, F1′)::Form1{X}
  F0::Form0{X}
  F0 == ⋆₀⁻¹{X}(∧₁₁′(F1, ⋆₁{X}(F1′)))
end
i1_imp′ = diag2dwd(I1Imp′, in_vars = [:F1, :F1′], out_vars = [:F0])
rules[:i₁′] = i1_imp′

max_stars_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_issues_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_forks_repo_name: AlgebraicJulia/DECAPODES-Benchmarks


text:
include("NonSimpleNumField/Field.jl")
include("NonSimpleNumField/Conjugates.jl")


text:
function solve()
    a, b = [parse(Int, x) for x in split(readline())]
    abs(a - b) < 3 ? "Yes" : "No"
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia