# 64. Optimal Taxation in an LQ Economy¶

Contents

## 64.1. Overview¶

In this lecture we study optimal fiscal policy in a linear quadratic setting.

We slightly modify a well-known model of Robert Lucas and Nancy Stokey [LS83] so that convenient formulas for solving linear-quadratic models can be applied to simplify the calculations.

The economy consists of a representative household and a benevolent government.

The government finances an exogenous stream of government purchases with state-contingent loans and a linear tax on labor income.

A linear tax is sometimes called a flat-rate tax.

The household maximizes utility by choosing paths for consumption and labor, taking prices and the government’s tax rate and borrowing plans as given.

Maximum attainable utility for the household depends on the government’s tax and borrowing plans.

The *Ramsey problem* [Ram27] is to choose tax and borrowing plans that maximize the household’s welfare, taking the household’s optimizing behavior as given.

There is a large number of competitive equilibria indexed by different government fiscal policies.

The Ramsey planner chooses the best competitive equilibrium.

We want to study the dynamics of tax rates, tax revenues, government debt under a Ramsey plan.

Because the Lucas and Stokey model features state-contingent government debt, the government debt dynamics differ substantially from those in a model of Robert Barro [Bar79].

The treatment given here closely follows this manuscript, prepared by Thomas J. Sargent and Francois R. Velde.

We cover only the key features of the problem in this lecture, leaving you to refer to that source for additional results and intuition.

### 64.1.1. Model Features¶

Linear quadratic (LQ) model

Representative household

Stochastic dynamic programming over an infinite horizon

Distortionary taxation

```
using LinearAlgebra, Statistics
```

## 64.2. The Ramsey Problem¶

We begin by outlining the key assumptions regarding technology, households and the government sector.

### 64.2.1. Technology¶

Labor can be converted one-for-one into a single, non-storable consumption good.

In the usual spirit of the LQ model, the amount of labor supplied in each period is unrestricted.

This is unrealistic, but helpful when it comes to solving the model.

Realistic labor supply can be induced by suitable parameter values.

### 64.2.2. Households¶

Consider a representative household who chooses a path \(\{\ell_t, c_t\}\) for labor and consumption to maximize

subject to the budget constraint

Here

\(\beta\) is a discount factor in \((0, 1)\)

\(p_t^0\) is a scaled Arrow-Debreu price at time \(0\) of history contingent goods at time \(t+j\)

\(b_t\) is a stochastic preference parameter

\(d_t\) is an endowment process

\(\tau_t\) is a flat tax rate on labor income

\(s_t\) is a promised time-\(t\) coupon payment on debt issued by the government

The scaled Arrow-Debreu price \(p^0_t\) is related to the unscaled Arrow-Debreu price as follows.

If we let \(\pi^0_t(x^t)\) denote the probability (density) of a history \(x^t = [x_t, x_{t-1}, \ldots, x_0]\) of the state \(x^t\), then the Arrow-Debreu time \(0\) price of a claim on one unit of consumption at date \(t\), history \(x^t\) would be

Thus, our scaled Arrow-Debreu price is the ordinary Arrow-Debreu price multiplied by the discount factor \(\beta^t\) and divided by an appropriate probability.

The budget constraint (64.2) requires that the present value of consumption be restricted to equal the present value of endowments, labor income and coupon payments on bond holdings.

### 64.2.3. Government¶

The government imposes a linear tax on labor income, fully committing to a stochastic path of tax rates at time zero.

The government also issues state-contingent debt.

Given government tax and borrowing plans, we can construct a competitive equilibrium with distorting government taxes.

Among all such competitive equilibria, the Ramsey plan is the one that maximizes the welfare of the representative consumer.

### 64.2.4. Exogenous Variables¶

Endowments, government expenditure, the preference shock process \(b_t\), and promised coupon payments on initial government debt \(s_t\) are all exogenous, and given by

\(d_t = S_d x_t\)

\(g_t = S_g x_t\)

\(b_t = S_b x_t\)

\(s_t = S_s x_t\)

The matrices \(S_d, S_g, S_b, S_s\) are primitives and \(\{x_t\}\) is an exogenous stochastic process taking values in \(\mathbb R^k\).

We consider two specifications for \(\{x_t\}\).

Discrete case: \(\{x_t\}\) is a discrete state Markov chain with transition matrix \(P\).

VAR case: \(\{x_t\}\) obeys \(x_{t+1} = A x_t + C w_{t+1}\) where \(\{w_t\}\) is independent zero mean Gaussian with identify covariance matrix.

### 64.2.5. Feasibility¶

The period-by-period feasibility restriction for this economy is

A labor-consumption process \(\{\ell_t, c_t\}\) is called *feasible* if (64.3) holds for all \(t\).

### 64.2.6. Government budget constraint¶

Where \(p_t^0\) is again a scaled Arrow-Debreu price, the time zero government budget constraint is

### 64.2.7. Equilibrium¶

An *equilibrium* is a feasible allocation \(\{\ell_t, c_t\}\), a sequence
of prices \(\{p_t^0\}\), and a tax system \(\{\tau_t\}\) such that

The allocation \(\{\ell_t, c_t\}\) is optimal for the household given \(\{p_t^0\}\) and \(\{\tau_t\}\).

The government’s budget constraint (64.4) is satisfied.

The *Ramsey problem* is to choose the equilibrium \(\{\ell_t, c_t, \tau_t, p_t^0\}\) that maximizes the
household’s welfare.

If \(\{\ell_t, c_t, \tau_t, p_t^0\}\) solves the Ramsey problem,
then \(\{\tau_t\}\) is called the *Ramsey plan*.

The solution procedure we adopt is

Use the first-order conditions from the household problem to pin down prices and allocations given \(\{\tau_t\}\).

Use these expressions to rewrite the government budget constraint (64.4) in terms of exogenous variables and allocations.

Maximize the household’s objective function (64.1) subject to the constraint constructed in step 2 and the feasibility constraint (64.3).

The solution to this maximization problem pins down all quantities of interest.

### 64.2.8. Solution¶

Step one is to obtain the first-conditions for the household’s problem, taking taxes and prices as given.

Letting \(\mu\) be the Lagrange multiplier on (64.2), the first-order conditions are \(p_t^0 = (c_t - b_t) / \mu\) and \(\ell_t = (c_t - b_t) (1 - \tau_t)\).

Rearranging and normalizing at \(\mu = b_0 - c_0\), we can write these conditions as

Substituting (64.5) into the government’s budget constraint (64.4) yields

The Ramsey problem now amounts to maximizing (64.1) subject to (64.6) and (64.3).

The associated Lagrangian is

The first order conditions associated with \(c_t\) and \(\ell_t\) are

and

Combining these last two equalities with (64.3) and working through the algebra, one can show that

where

\(\nu := \lambda / (1 + 2 \lambda)\)

\(\bar \ell_t := (b_t - d_t + g_t) / 2\)

\(\bar c_t := (b_t + d_t - g_t) / 2\)

\(m_t := (b_t - d_t - s_t ) / 2\)

Apart from \(\nu\), all of these quantities are expressed in terms of exogenous variables.

To solve for \(\nu\), we can use the government’s budget constraint again.

The term inside the brackets in (64.6) is \((b_t - c_t)(s_t + g_t) - (b_t - c_t) \ell_t + \ell_t^2\).

Using (64.8), the definitions above and the fact that \(\bar \ell = b - \bar c\), this term can be rewritten as

Reinserting into (64.6), we get

Although it might not be clear yet, we are nearly there because:

The two expectations terms in (64.9) can be solved for in terms of model primitives.

This in turn allows us to solve for the Lagrange multiplier \(\nu\).

With \(\nu\) in hand, we can go back and solve for the allocations via (64.8).

Once we have the allocations, prices and the tax system can be derived from (64.5).

### 64.2.9. Computing the Quadratic Term¶

Let’s consider how to obtain the term \(\nu\) in (64.9).

If we can compute the two expected geometric sums

then the problem reduces to solving

for \(\nu\).

Provided that \(4 b_0 < a_0\), there is a unique solution \(\nu \in (0, 1/2)\), and a unique corresponding \(\lambda > 0\).

Let’s work out how to compute mathematical expectations in (64.10).

For the first one, the random variable \((b_t - \bar c_t) (g_t + s_t )\) inside the summation can be expressed as

For the second expectation in (64.10), the random variable \(2 m_t^2\) can be written as

It follows that both objects of interest are special cases of the expression

where \(H\) is a matrix conformable to \(x_t\) and \(x_t'\) is the transpose of column vector \(x_t\).

Suppose first that \(\{x_t\}\) is the Gaussian VAR described above.

In this case, the formula for computing \(q(x_0)\) is known to be \(q(x_0) = x_0' Q x_0 + v\), where

\(Q\) is the solution to \(Q = H + \beta A' Q A\), and

\(v = \text{trace} \, (C' Q C) \beta / (1 - \beta)\)

The first equation is known as a discrete Lyapunov equation, and can be solved using this function.

### 64.2.10. Finite state Markov case¶

Next suppose that \(\{x_t\}\) is the discrete Markov process described above.

Suppose further that each \(x_t\) takes values in the state space \(\{x^1, \ldots, x^N\} \subset \mathbb R^k\).

Let \(h \colon \mathbb R^k \to \mathbb R\) be a given function, and suppose that we wish to evaluate

For example, in the discussion above, \(h(x_t) = x_t' H x_t\).

It is legitimate to pass the expectation through the sum, leading to

Here

\(P^t\) is the \(t\)-th power of the transition matrix \(P\)

\(h\) is, with some abuse of notation, the vector \((h(x^1), \ldots, h(x^N))\)

\((P^t h)[j]\) indicates the \(j\)-th element of \(P^t h\)

It can be show that (64.12) is in fact equal to the \(j\)-th element of the vector \((I - \beta P)^{-1} h\).

This last fact is applied in the calculations below.

### 64.2.11. Other Variables¶

We are interested in tracking several other variables besides the ones described above.

To prepare the way for this, we define

as the scaled Arrow-Debreu time \(t\) price of a history contingent claim on one unit of consumption at time \(t+j\).

These are prices that would prevail at time \(t\) if market were reopened at time \(t\).

These prices are constituents of the present value of government obligations outstanding at time \(t\), which can be expressed as

Using our expression for prices and the Ramsey plan, we can also write \(B_t\) as

This version is more convenient for computation.

Using the equation

it is possible to verity that (64.13) implies that

and

Define

\(R_{t}\) is the gross \(1\)-period risk-free rate for loans between \(t\) and \(t+1\).

### 64.2.12. A Martingale¶

We now want to study the following two objects, namely,

and the cumulation of \(\pi_t\)

The term \(\pi_{t+1}\) is the difference between two quantities:

\(B_{t+1}\), the value of government debt at the start of period \(t+1\).

\(R_t [B_t + g_t - \tau_t ]\), which is what the government would have owed at the beginning of period \(t+1\) if it had simply borrowed at the one-period risk-free rate rather than selling state-contingent securities.

Thus, \(\pi_{t+1}\) is the excess payout on the actual portfolio of state contingent government debt relative to an alternative portfolio sufficient to finance \(B_t + g_t - \tau_t \ell_t\) and consisting entirely of risk-free one-period bonds.

Use expressions (64.14) and (64.15) to obtain

or

where \(\tilde E_t\) is the conditional mathematical expectation taken with respect to a one-step transition density that has been formed by multiplying the original transition density with the likelihood ratio

It follows from equation (64.16) that

which asserts that \(\{\pi_{t+1}\}\) is a martingale difference sequence under the distorted probability measure, and that \(\{\Pi_t\}\) is a martingale under the distorted probability measure.

In the tax-smoothing model of Robert Barro [Bar79], government debt is a random walk.

In the current model, government debt \(\{B_t\}\) is not a random walk, but the `excess payoff`

\(\{\Pi_t\}\) on it is.

## 64.3. Implementation¶

The following code provides functions for

Solving for the Ramsey plan given a specification of the economy.

Simulating the dynamics of the major variables.

Description and clarifications are given below

```
using LaTeXStrings, QuantEcon, Plots, LinearAlgebra, Parameters
abstract type AbstractStochProcess end
struct ContStochProcess{TF <: AbstractFloat} <: AbstractStochProcess
A::Matrix{TF}
C::Matrix{TF}
end
struct DiscreteStochProcess{TF <: AbstractFloat} <: AbstractStochProcess
P::Matrix{TF}
x_vals::Matrix{TF}
end
struct Economy{TF <: AbstractFloat, SP <: AbstractStochProcess}
β::TF
Sg::Matrix{TF}
Sd::Matrix{TF}
Sb::Matrix{TF}
Ss::Matrix{TF}
proc::SP
end
function compute_exog_sequences(econ, x)
# compute exogenous variable sequences
Sg, Sd, Sb, Ss = econ.Sg, econ.Sd, econ.Sb, econ.Ss
g, d, b, s = [dropdims(S * x, dims = 1) for S in (Sg, Sd, Sb, Ss)]
#= solve for Lagrange multiplier in the govt budget constraint
In fact we solve for ν = λ / (1 + 2*λ). Here ν is the
solution to a quadratic equation a(ν^2 - ν) + b = 0 where
a and b are expected discounted sums of quadratic forms of the state. =#
Sm = Sb - Sd - Ss
return g, d, b, s, Sm
end
function compute_allocation(econ, Sm, ν, x, b)
Sg, Sd, Sb, Ss = econ.Sg, econ.Sd, econ.Sb, econ.Ss
# solve for the allocation given ν and x
Sc = 0.5 .* (Sb + Sd - Sg - ν .* Sm)
Sl = 0.5 .* (Sb - Sd + Sg - ν .* Sm)
c = dropdims(Sc * x, dims = 1)
l = dropdims(Sl * x, dims = 1)
p = dropdims((Sb - Sc) * x, dims = 1) # Price without normalization
τ = 1 .- l ./ (b .- c)
rvn = l .* τ
return Sc, Sl, c, l, p, τ, rvn
end
function compute_ν(a0, b0)
disc = a0^2 - 4a0 * b0
if disc ≥ 0
ν = 0.5 *(a0 - sqrt(disc)) / a0
else
println("There is no Ramsey equilibrium for these parameters.")
error("Government spending (economy.g) too low")
end
# Test that the Lagrange multiplier has the right sign
if ν * (0.5 - ν) < 0
print("Negative multiplier on the government budget constraint.")
error("Government spending (economy.g) too low")
end
return ν
end
function compute_Π(B, R, rvn, g, ξ)
π = B[2:end] - R[1:end-1] .* B[1:end-1] - rvn[1:end-1] + g[1:end-1]
Π = cumsum(π .* ξ)
return π, Π
end
function compute_paths(econ::Economy{<:AbstractFloat, <:DiscreteStochProcess}, T)
# simplify notation
@unpack β, Sg, Sd, Sb, Ss = econ
@unpack P, x_vals = econ.proc
mc = MarkovChain(P)
state = simulate(mc, T, init=1)
x = x_vals[:, state]
# Compute exogenous sequence
g, d, b, s, Sm = compute_exog_sequences(econ, x)
# compute a0, b0
ns = size(P, 1)
F = I - β.*P
a0 = (F \ ((Sm * x_vals)'.^2))[1] ./ 2
H = ((Sb - Sd + Sg) * x_vals) .* ((Sg - Ss)*x_vals)
b0 = (F \ H')[1] ./ 2
# compute lagrange multiplier
ν = compute_ν(a0, b0)
# Solve for the allocation given ν and x
Sc, Sl, c, l, p, τ, rvn = compute_allocation(econ, Sm, ν, x, b)
# compute remaining variables
H = ((Sb - Sc) * x_vals) .* ((Sl - Sg) * x_vals) - (Sl * x_vals).^2
temp = dropdims(F * H', dims = 2)
B = temp[state] ./ p
H = dropdims(P[state, :] * ((Sb - Sc) * x_vals)', dims = 2)
R = p ./ (β .* H)
temp = dropdims(P[state, :] *((Sb - Sc) * x_vals)', dims = 2)
ξ = p[2:end] ./ temp[1:end-1]
# compute π
π, Π = compute_Π(B, R, rvn, g, ξ)
return (g = g, d = d, b = b, s = s, c = c,
l = l, p = p, τ = τ, rvn = rvn, B = B,
R = R, π = π, Π = Π, ξ = ξ)
end
function compute_paths(econ::Economy{<:AbstractFloat, <:ContStochProcess}, T)
# simplify notation
@unpack β, Sg, Sd, Sb, Ss = econ
@unpack A, C = econ.proc
# generate an initial condition x0 satisfying x0 = A x0
nx, nx = size(A)
x0 = nullspace(I - A)
x0 = x0[end] < 0 ? -x0 : x0
x0 = x0 ./ x0[end]
x0 = dropdims(x0, dims = 2)
# generate a time series x of length T starting from x0
nx, nw = size(C)
x = zeros(nx, T)
w = randn(nw, T)
x[:, 1] = x0
for t in 2:T
x[:, t] = A *x[:, t-1] + C * w[:, t]
end
# compute exogenous sequence
g, d, b, s, Sm = compute_exog_sequences(econ, x)
# compute a0 and b0
H = Sm'Sm
a0 = 0.5 * var_quadratic_sum(A, C, H, β, x0)
H = (Sb - Sd + Sg)'*(Sg + Ss)
b0 = 0.5 * var_quadratic_sum(A, C, H, β, x0)
# compute lagrange multiplier
ν = compute_ν(a0, b0)
# solve for the allocation given ν and x
Sc, Sl, c, l, p, τ, rvn = compute_allocation(econ, Sm, ν, x, b)
# compute remaining variables
H = Sl'Sl - (Sb - Sc)' *(Sl - Sg)
L = zeros(T)
for t in eachindex(L)
L[t] = var_quadratic_sum(A, C, H, β, x[:, t])
end
B = L ./ p
Rinv = dropdims(β .* (Sb- Sc)*A*x, dims = 1) ./ p
R = 1 ./ Rinv
AF1 = (Sb - Sc) * x[:, 2:end]
AF2 = (Sb - Sc) * A * x[:, 1:end-1]
ξ = AF1 ./ AF2
ξ = dropdims(ξ, dims = 1)
# compute π
π, Π = compute_Π(B, R, rvn, g, ξ)
return (g = g, d = d, b = b, s = s,
c = c, l = l, p = p, τ = τ,
rvn = rvn, B = B, R = R,
π = π, Π = Π, ξ = ξ)
end
function gen_fig_1(path)
T = length(path.c)
plt_1 = plot(path.rvn, lw=2, label = L"\tau_t l_t")
plot!(plt_1, path.g, lw=2, label= L"g_t")
plot!(plt_1, path.c, lw=2, label= L"c_t")
plot!(xlabel="Time", grid=true)
plt_2 = plot(path.rvn, lw=2, label=L"\tau_t l_t")
plot!(plt_2, path.g, lw=2, label=L"g_t")
plot!(plt_2, path.B[2:end], lw=2, label=L"B_{t+1}")
plot!(xlabel="Time", grid=true)
plt_3 = plot(path.R, lw=2, label=L"R_{t-1}")
plot!(plt_3, xlabel="Time", grid=true)
plt_4 = plot(path.rvn, lw=2, label=L"\tau_t l_t")
plot!(plt_4, path.g, lw=2, label=L"g_t")
plot!(plt_4, path.π, lw=2, label=L"\pi_t")
plot!(plt_4, xlabel="Time", grid=true)
plot(plt_1, plt_2, plt_3, plt_4, layout=(2,2), size = (800,600))
end
function gen_fig_2(path)
T = length(path.c)
paths = [path.ξ, path.Π]
labels = [L"\xi_t", L"\Pi_t"]
plt_1 = plot()
plt_2 = plot()
plots = [plt_1, plt_2]
for (plot, path, label) in zip(plots, paths, labels)
plot!(plot, 2:T, path, lw=2, label=label, xlabel="Time", grid=true)
end
plot(plt_1, plt_2, layout=(2,1), size = (600,500))
end
```

```
gen_fig_2 (generic function with 1 method)
```

### 64.3.1. Comments on the Code¶

The function `var_quadratic_sum`

From `QuantEcon.jl`

is for computing the value of (64.11)
when the exogenous process \(\{ x_t \}\) is of the VAR type described above.

This code defines two Types: `Economy`

and `Path`

.

The first is used to collect all the parameters and primitives of a given LQ economy, while the second collects output of the computations.

## 64.4. Examples¶

Let’s look at two examples of usage.

### 64.4.1. The Continuous Case¶

Our first example adopts the VAR specification described above.

Regarding the primitives, we set

\(\beta = 1 / 1.05\)

\(b_t = 2.135\) and \(s_t = d_t = 0\) for all \(t\)

Government spending evolves according to

with \(\rho = 0.7\), \(\mu_g = 0.35\) and \(C_g = \mu_g \sqrt{1 - \rho^2} / 10\).

Here’s the code

```
# for reproducible results
using Random
Random.seed!(42)
# parameters
β = 1 / 1.05
ρ, mg = .7, .35
A = [ρ mg*(1 - ρ); 0.0 1.0]
C = [sqrt(1 - ρ^2) * mg / 10 0.0; 0 0]
Sg = [1.0 0.0]
Sd = [0.0 0.0]
Sb = [0 2.135]
Ss = [0.0 0.0]
proc = ContStochProcess(A, C)
econ = Economy(β, Sg, Sd, Sb, Ss, proc)
T = 50
path = compute_paths(econ, T)
gen_fig_1(path)
```

The legends on the figures indicate the variables being tracked.

Most obvious from the figure is tax smoothing in the sense that tax revenue is much less variable than government expenditure

```
gen_fig_2(path)
```

See the original manuscript for comments and interpretation

### 64.4.2. The Discrete Case¶

Our second example adopts a discrete Markov specification for the exogenous process

```
# Parameters
β = 1 / 1.05
P = [0.8 0.2 0.0
0.0 0.5 0.5
0.0 0.0 1.0]
# Possible states of the world
# Each column is a state of the world. The rows are [g d b s 1]
x_vals = [0.5 0.5 0.25;
0.0 0.0 0.0;
2.2 2.2 2.2;
0.0 0.0 0.0;
1.0 1.0 1.0]
Sg = [1.0 0.0 0.0 0.0 0.0]
Sd = [0.0 1.0 0.0 0.0 0.0]
Sb = [0.0 0.0 1.0 0.0 0.0]
Ss = [0.0 0.0 0.0 1.0 0.0]
proc = DiscreteStochProcess(P, x_vals)
econ = Economy(β, Sg, Sd, Sb, Ss, proc)
T = 15
path = compute_paths(econ, T)
gen_fig_1(path)
```

The call `gen_fig_2(path)`

generates

```
gen_fig_2(path)
```

See the original manuscript for comments and interpretation

## 64.5. Exercises¶

### 64.5.1. Exercise 1¶

Modify the VAR example given above, setting

with \(\rho = 0.95\) and \(C_g = 0.7 \sqrt{1 - \rho^2}\).

Produce the corresponding figures.

## 64.6. Solutions¶

```
# parameters
β = 1 / 1.05
ρ, mg = .95, .35
A = [0. 0. 0. ρ mg*(1-ρ);
1. 0. 0. 0. 0.;
0. 1. 0. 0. 0.;
0. 0. 1. 0. 0.;
0. 0. 0. 0. 1.]
C = zeros(5, 5)
C[1, 1] = sqrt(1 - ρ^2) * mg / 8
Sg = [1. 0. 0. 0. 0.]
Sd = [0. 0. 0. 0. 0.]
Sb = [0. 0. 0. 0. 2.135]
Ss = [0. 0. 0. 0. 0.]
proc = ContStochProcess(A, C)
econ = Economy(β, Sg, Sd, Sb, Ss, proc)
T = 50
path = compute_paths(econ, T)
```

```
(g = [0.3500000000000022, 0.35651269112333617, 0.33872515990993246, 0.3697704287947342, 0.3311095154846114, 0.3650509399142057, 0.34162861887273244, 0.3627748524207786, 0.33538802774966137, 0.35316315802093146 … 0.3868534729995743, 0.35428035513450773, 0.3265551424807309, 0.3897782718961585, 0.38813503098337354, 0.33290449602825256, 0.3389345173160616, 0.3741739007644913, 0.38469878165606924, 0.3070933401453645], d = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], b = [2.135, 2.135, 2.135, 2.135, 2.135, 2.135, 2.135, 2.135, 2.135, 2.135 … 2.135, 2.135, 2.135, 2.135, 2.135, 2.135, 2.135, 2.135, 2.135, 2.135], s = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], c = [0.6181319813258108, 0.6148756357641438, 0.6237694013708457, 0.6082467669284448, 0.6275772235835062, 0.610606511368709, 0.6223176718894456, 0.6117445551154226, 0.6254379674509811, 0.6165504023153461 … 0.5997052448260247, 0.615991803758558, 0.6298544100854464, 0.5982428453777326, 0.5990644658341251, 0.6266797333116856, 0.6236647226677811, 0.6060450309435662, 0.6007825904977773, 0.6395853112531296], l = [0.9681319813258129, 0.9713883268874799, 0.962494561280778, 0.9780171957231789, 0.9586867390681175, 0.9756574512829147, 0.9639462907621781, 0.9745194075362011, 0.9608259952006426, 0.9697135603362776 … 0.986558717825599, 0.9702721588930657, 0.9564095525661773, 0.9880211172738911, 0.9871994968174986, 0.9595842293399381, 0.9625992399838427, 0.9802189317080575, 0.9854813721538465, 0.9466786513984942], p = [1.516868018674189, 1.520124364235856, 1.511230598629154, 1.526753233071555, 1.5074227764164936, 1.5243934886312906, 1.5126823281105541, 1.5232554448845772, 1.5095620325490184, 1.5184495976846535 … 1.5352947551739748, 1.5190081962414417, 1.5051455899145534, 1.5367571546222671, 1.5359355341658745, 1.508320266688314, 1.5113352773322186, 1.5289549690564335, 1.5342174095022225, 1.49541468874687], τ = [0.36175595410601125, 0.36098101593432297, 0.36310543066434575, 0.3594137058052381, 0.36402265239275045, 0.35997007428906724, 0.36275695640193384, 0.3602390125642754, 0.36350678244191825, 0.3613791581789043 … 0.35741412878479806, 0.36124626496824774, 0.36457339477673223, 0.3570740085366674, 0.35726501870821026, 0.3638060493300852, 0.3630802811120739, 0.35889614047104235, 0.3576651092275206, 0.36694573182787626], rvn = [0.3502275086050625, 0.35065274510658473, 0.3494870021859474, 0.35151278465611463, 0.3489836895693328, 0.3512074852189928, 0.3496782225718213, 0.3510599090955638, 0.34926676600193957, 0.35043427010899214 … 0.3526100246266839, 0.3505071934027982, 0.3486814773759468, 0.352796660863865, 0.3526918466992394, 0.3491025474756173, 0.34950080265160227, 0.3517967914066701, 0.35247230261309237, 0.3473796905432474], B = [2.121418612773956e-15, 0.02785687092981279, -0.025259330707965846, 0.06922571840826187, -0.02690101331991359, 0.028363380073783643, 0.025104498102512453, 0.008499238954541397, 0.01054561187869878, -0.022580577921653536 … 0.19556886178896726, 0.11338120517671109, 0.15618372453125268, 0.24714211635220124, 0.23403960104186594, 0.12141181447963512, 0.19160384679514283, 0.15073523800547306, 0.12747081968262972, -0.011310843257403065], R = [1.0500000000000003, 1.052254093828637, 1.046097688807191, 1.0568427015333257, 1.0413381301686913, 1.0589480372948412, 1.0406598573954826, 1.0606959876256963, 1.0400408331592836, 1.0538574371418683 … 1.056600300617863, 1.0630260588571803, 1.0334005156376416, 1.051631234650584, 1.0617756601230834, 1.051805091227126, 1.0332989878230514, 1.0458771189834337, 1.067725457865234, 1.0387490388418728], π = [0.027629362324750217, -0.048711891168365884, 0.08488760360666647, -0.08180406443946675, 0.038502256899264176, 0.008912607137544493, -0.025675608159893226, 0.013245446546959527, -0.04742718313842803, -0.03661141937371698 … -0.03624855295757973, -0.05901346460811463, 0.03942971057548572, 0.06361543999218144, 0.011118843120535737, -0.09164255312740105, 0.047704230742974585, -0.05781510828649172, -0.0078026074139954105, -0.1151882035245187], Π = [0.027688675824850597, -0.020842178240535077, 0.06459862587599774, -0.016696060557747566, 0.02191846217956882, 0.030837967367404434, 0.005212886277779074, 0.018472975806920425, -0.028780943926234998, -0.06528992928746707 … 0.2470192240760901, 0.18826475664273654, 0.2278193199622608, 0.29174401184596843, 0.3028741748450268, 0.21187001912829437, 0.25975178070392846, 0.2021929552249885, 0.19439423517562407, 0.08022396221353403], ξ = [1.0021467560272734, 0.9962835131497055, 1.0065168586031672, 0.9937731944101827, 1.002915824865694, 1.0007739654833503, 0.9980320984043216, 1.0011055106469917, 0.996346748977967, 0.9972021294383785 … 0.9973319838035288, 0.9956112189568775, 1.0031664636188369, 1.0048612709676181, 1.0010180806042441, 0.9930338321131117, 1.0037214903142673, 0.9955672000771528, 0.9995017864638446, 0.9911628922816564])
```

```
gen_fig_1(path)
```

```
gen_fig_2(path)
```