8. Markov Jump Linear Quadratic Dynamic Programming#

In addition to what’s in Anaconda, this lecture will need the following libraries:

!pip install --upgrade quantecon

8.1. Overview#

This lecture describes Markov jump linear quadratic dynamic programming, an extension of the method described in the first LQ control lecture.

Markov jump linear quadratic dynamic programming is described and analyzed in [Do Val et al., 1999] and the references cited there.

The method has been applied to problems in macroeconomics and monetary economics by [Svensson et al., 2008] and [Svensson and Williams, 2009].

The periodic models of seasonality described in chapter 14 of [Hansen and Sargent, 2013] are a special case of Markov jump linear quadratic problems.

Markov jump linear quadratic dynamic programming combines advantages of

the computational simplicity of linear quadratic dynamic programming, with
the ability of finite state Markov chains to represent interesting patterns of random variation.

The idea is to replace the constant matrices that define a linear quadratic dynamic programming problem with \(N\) sets of matrices that are fixed functions of the state of an \(N\) state Markov chain.

The state of the Markov chain together with the continuous \(n \times 1\) state vector \(x_t\) form the state of the system.

For the class of infinite horizon problems being studied in this lecture, we obtain \(N\) interrelated matrix Riccati equations that determine \(N\) optimal value functions and \(N\) linear decision rules.

One of these value functions and one of these decision rules apply in each of the \(N\) Markov states.

That is, when the Markov state is in state \(j\), the value function and the decision rule for state \(j\) prevails.

8.2. Review of useful LQ dynamic programming formulas#

To begin, it is handy to have the following reminder in mind.

A linear quadratic dynamic programming problem consists of a scalar discount factor \(\beta \in (0,1)\), an \(n\times 1\) state vector \(x_t\), an initial condition for \(x_0\), a \(k \times 1\) control vector \(u_t\), a \(p \times 1\) random shock vector \(w_{t+1}\) and the following two triples of matrices:

A triple of matrices \((R, Q, W)\) defining a loss function

\[ r(x_t, u_t) = x_t' R x_t + u_t' Q u_t + 2 u_t' W x_t \]

a triple of matrices \((A, B, C)\) defining a state-transition law

\[ x_{t+1} = A x_t + B u_t + C w_{t+1} \]

The problem is

\[ -x_0' P x_0 - \rho = \min_{\{u_t\}_{t=0}^\infty} E \sum_{t=0}^{\infty} \beta^t r(x_t, u_t) \]

subject to the transition law for the state.

The optimal decision rule has the form

\[ u_t = - F x_t \]

and the optimal value function is of the form

\[ -\left( x_t' P x_t + \rho \right) \]

where \(P\) solves the algebraic matrix Riccati equation

\[ P = R+ \beta A' P A -(\beta B' P A + W)' (Q + \beta B P B )^{-1} (\beta B P A + W) \]

and the constant \(\rho\) satisfies

\[ \rho = \beta \left( \rho + {\rm trace}(P C C') \right) \]

and the matrix \(F\) in the decision rule for \(u_t\) satisfies

\[ F = (Q + \beta B' P B)^{-1} (\beta (B' P A )+ W) \]

With the preceding formulas in mind, we are ready to approach Markov Jump linear quadratic dynamic programming.

8.3. Linked Riccati equations for Markov LQ dynamic programming#

The key idea is to make the matrices \(A, B, C, R, Q, W\) fixed functions of a finite state \(s\) that is governed by an \(N\) state Markov chain.

This makes decision rules depend on the Markov state, and so fluctuate through time in limited ways.

In particular, we use the following extension of a discrete-time linear quadratic dynamic programming problem.

We let \(s_t \in [1, 2, \ldots, N]\) be a time \(t\) realization of an \(N\)-state Markov chain with transition matrix \(\Pi\) having typical element \(\Pi_{ij}\).

Here \(i\) denotes today and \(j\) denotes tomorrow and

\[ \Pi_{ij} = {\rm Prob}(s_{t+1} = j |s_t = i) \]

We’ll switch between labeling today’s state as \(s_t\) and \(i\) and between labeling tomorrow’s state as \(s_{t+1}\) or \(j\).

The decision-maker solves the minimization problem:

\[ \min_{\{u_t\}_{t=0}^\infty} E \sum_{t=0}^{\infty} \beta^t r(x_t, s_t, u_t) \]

with

\[ r(x_t, s_t, u_t) = x_t' R_{s_t} x_t + u_t' Q_{s_t} u_t + 2 u_t' W_{s_t} x_t \]

subject to linear laws of motion with matrices \((A,B,C)\) each possibly dependent on the Markov-state-\(s_t\):

\[ x_{t+1} = A_{s_t} x_t + B_{s_t} u_t + C_{s_t} w_{t+1} \]

where \(\{w_{t+1}\}\) is an i.i.d. stochastic process with \(w_{t+1} \sim {\cal N}(0,I)\).

The optimal decision rule for this problem has the form

\[ u_t = - F_{s_t} x_t \]

and the optimal value functions are of the form

\[ -\left( x_t' P_{s_t} x_t + \rho_{s_t} \right) \]

or equivalently

\[ -x_t' P_i x_t - \rho_i \]

The optimal value functions \(- x' P_i x - \rho_i\) for \(i = 1, \ldots, n\) satisfy the \(N\) interrelated Bellman equations

\[\begin{split} \begin{split} -x' P_i x - \rho_i & = \max_u - \\ & \left[ x'R_i x + u' Q_i u + 2 u' W_i x + \beta \sum_j \Pi_{ij}E ((A_i x + B_i u + C_i w)' P_j (A_i x + B_i u + C_i w) x + \rho_j) \right] \end{split} \end{split}\]

The matrices \(P_{s_t} = P_i\) and the scalars \(\rho_{s_t} = \rho_i, i = 1, \ldots\), n satisfy the following stacked system of algebraic matrix Riccati equations:

\[ P_i = R_i + \beta \sum_j A_i' P_j A_i \Pi_{ij} -\sum_j \Pi_{ij}[ (\beta B_i' P_j A_i + W_i)' (Q + \beta B_i' P_j B_i)^{-1} (\beta B_i' P_j A_i + W_i)] \]

\[ \rho_i = \beta \sum_j \Pi_{ij} ( \rho_j + {\rm trace}(P_j C_i C_i') ) \]

and the \(F_i\) in the optimal decision rules are

\[ F_i = (Q_i + \beta \sum_j \Pi_{ij} B_i' P_j B_i)^{-1} (\beta \sum_j \Pi_{ij}(B_i' P_j A_i )+ W_i) \]

8.4. Applications#

We now describe some Python code and a few examples that put the code to work.

To begin, we import these Python modules

import numpy as np
import quantecon as qe
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Set discount factor
β = 0.95

8.5. Example 1#

This example is a version of a classic problem of optimally adjusting a variable \(k_t\) to a target level in the face of costly adjustment.

This provides a model of gradual adjustment.

Given \(k_{0}\), the objective function is

\[ \max_{\left\{ k_{t}\right\} _{t=1}^{\infty}}E_{0}\sum_{t=0}^{\infty}\beta^{t}r\left(s_{t},k_{t}\right) \]

where the one-period payoff function is

\[ r(s_{t},k_{t})=f_{1,s_t} k_{t}-f_{2,s_t} k_{t}^{2}-d_{s_t}(k_{t+1}-k_{t})^{2}, \]

\(E_0\) is a mathematical expectation conditioned on time \(0\) information \(x_0, s_0\) and the transition law for continuous state variable \(k_t\) is

\[ k_{t+1}-k_{t}=u_{t} \]

We can think of \(k_t\) as the decision-maker’s capital and \(u_t\) as costs of adjusting the level of capital.

We assume that \(f_{1}\left(s_{t}\right)>0\), \(f_{2}\left(s_{t}\right)>0\), and \(d\left(s_{t}\right)>0\).

Denote the state transition matrix for Markov state \(s_{t}\in\left\{1,2 \right\}\) as \(\Pi\):

\[ \Pr \left(s_{t+1}=j \mid s_{t}= i \right)=\Pi_{ij} \]

Let \(x_{t}=\begin{bmatrix} k_{t}\\ 1 \end{bmatrix}\)

We can represent the one-period payoff function \(r\left(s_{t},k_{t}\right)\) and the state-transition law as

\[\begin{split} \begin{aligned} r\left(s_{t},k_{t}\right) =f_{1,s_t} k_{t}-f_{2,s_t} k_{t}^{2}-d_{s_t} u_{t}{}^{2} \\ =-x_{t}^{\prime}\underset{\equiv R(s_{t})}{\underbrace{\begin{bmatrix} f_{2,s_t} & -\frac{f_{1,s_t}}{2}\\ -\frac{f_{1,s_t}}{2} & 0 \end{bmatrix}}}x_{t}+\underset{\equiv Q\left(s_{t}\right)}{\underbrace{d_{s_t}}}u_{t}{}^{2} \end{aligned} \end{split}\]

\[\begin{split} x_{t+1}=\begin{bmatrix} k_{t+1}\\ 1 \end{bmatrix}=\underset{\equiv A\left(s_{t}\right)}{\underbrace{I_{2}}}x_{t}+\underset{\equiv B\left(s_{t}\right)}{\underbrace{\begin{bmatrix} 1\\ 0 \end{bmatrix}}}u_{t} \end{split}\]

def construct_arrays1(f1_vals=[1. ,1.],
                      f2_vals=[1., 1.],
                      d_vals=[1., 1.]):
    """
    Construct matrices that map the problem described in example 1
    into a Markov jump linear quadratic dynamic programming problem
    """

    # Number of Markov states
    m = len(f1_vals)
    # Number of state and control variables
    n, k = 2, 1

    # Construct sets of matrices for each state
    As = [np.eye(n) for i in range(m)]
    Bs = [np.array([[1, 0]]).T for i in range(m)]

    Rs = np.zeros((m, n, n))
    Qs = np.zeros((m, k, k))

    for i in range(m):
        Rs[i, 0, 0] = f2_vals[i]
        Rs[i, 1, 0] = - f1_vals[i] / 2
        Rs[i, 0, 1] = - f1_vals[i] / 2

        Qs[i, 0, 0] = d_vals[i]

    Cs, Ns = None, None

    # Compute the optimal k level of the payoff function in each state
    k_star = np.empty(m)
    for i in range(m):
        k_star[i] = f1_vals[i] / (2 * f2_vals[i])

    return Qs, Rs, Ns, As, Bs, Cs, k_star

The continuous part of the state \(x_t\) consists of two variables, namely, \(k_t\) and a constant term.

state_vec1 = ["k", "constant term"]

We start with a Markov transition matrix that makes the Markov state be strictly periodic:

\[\begin{split} \Pi_{1}=\begin{bmatrix} 0 & 1\\ 1 & 0 \end{bmatrix}, \end{split}\]

We set \(f_{1,{s_t}}\) and \(f_{2,{s_t}}\) to be independent of the Markov state \(s_t\)

\[ f_{1,1}=f_{1,2} = 1, \]

\[ f_{2,1} =f_{2,2} = 1 \]

In contrast to \(f_{1,{s_t}}\) and \(f_{2,{s_t}}\), we make the adjustment cost \(d_{s_t}\) vary across Markov states \(s_t\).

We set the adjustment cost to be lower in Markov state \(2\)

\[ d_1=1, d_2 = 0.5 \]

The following code forms a Markov switching LQ problem and computes the optimal value functions and optimal decision rules for each Markov state

# Construct Markov transition matrix
Π1 = np.array([[0., 1.],
               [1., 0.]])

# Construct matrices
Qs, Rs, Ns, As, Bs, Cs, k_star = construct_arrays1(d_vals=[1., 0.5])

# Construct a Markov Jump LQ problem
ex1_a = qe.LQMarkov(Π1, Qs, Rs, As, Bs, Cs=Cs, Ns=Ns, beta=β)
# Solve for optimal value functions and decision rules
ex1_a.stationary_values();

Let’s look at the value function matrices and the decision rules for each Markov state

# P(s)
ex1_a.Ps

array([[[ 1.56626026, -0.78313013],
        [-0.78313013, -4.60843493]],

       [[ 1.37424214, -0.68712107],
        [-0.68712107, -4.65643947]]])

# d(s) = 0, since there is no randomness
ex1_a.ds

array([0., 0.])

# F(s)
ex1_a.Fs

array([[[ 0.56626026, -0.28313013]],

       [[ 0.74848427, -0.37424214]]])

Now we’ll plot the decision rules and see if they make sense

# Plot the optimal decision rules
k_grid = np.linspace(0., 1., 100)
# Optimal choice in state s1
u1_star = - ex1_a.Fs[0, 0, 1] - ex1_a.Fs[0, 0, 0] * k_grid
# Optimal choice in state s2
u2_star = - ex1_a.Fs[1, 0, 1] - ex1_a.Fs[1, 0, 0] * k_grid

fig, ax = plt.subplots()
ax.plot(k_grid, k_grid + u1_star, label="$\overline{s}_1$ (high)")
ax.plot(k_grid, k_grid + u2_star, label="$\overline{s}_2$ (low)")

# The optimal k*
ax.scatter([0.5, 0.5], [0.5, 0.5], marker="*")
ax.plot([k_star[0], k_star[0]], [0., 1.0], '--')

# 45 degree line
ax.plot([0., 1.], [0., 1.], '--', label="45 degree line")

ax.set_xlabel("$k_t$")
ax.set_ylabel("$k_{t+1}$")
ax.legend()
plt.show()

<>:9: SyntaxWarning: invalid escape sequence '\o'
<>:10: SyntaxWarning: invalid escape sequence '\o'
<>:9: SyntaxWarning: invalid escape sequence '\o'
<>:10: SyntaxWarning: invalid escape sequence '\o'
/tmp/ipykernel_7172/454231666.py:9: SyntaxWarning: invalid escape sequence '\o'
  ax.plot(k_grid, k_grid + u1_star, label="$\overline{s}_1$ (high)")
/tmp/ipykernel_7172/454231666.py:10: SyntaxWarning: invalid escape sequence '\o'
  ax.plot(k_grid, k_grid + u2_star, label="$\overline{s}_2$ (low)")

_images/361fde00e13229ed65fd01d8a70391dcd2804d8ddd5f3832edc2ba3014fa66a5.png

The above graph plots \(k_{t+1}= k_t + u_t = k_t - F x_t\) as an affine (i.e., linear in \(k_t\) plus a constant) function of \(k_t\) for both Markov states \(s_t\).

It also plots the 45 degree line.

Notice that the two \(s_t\)-dependent closed loop functions that determine \(k_{t+1}\) as functions of \(k_t\) share the same rest point (also called a fixed point) at \(k_t = 0.5\).

Evidently, the optimal decision rule in Markov state \(2\), in which the adjustment cost is lower, makes \(k_{t+1}\) a flatter function of \(k_t\) in Markov state \(2\).

This happens because when \(k_t\) is not at its fixed point, \(| u_{t,2}| > | u_{t,2} |\), so that the decision-maker adjusts toward the fixed point faster when the Markov state \(s_t\) takes a value that makes it cheaper.

# Compute time series
T = 20
x0 = np.array([[0., 1.]]).T
x_path = ex1_a.compute_sequence(x0, ts_length=T)[0]

fig, ax = plt.subplots()
ax.plot(range(T), x_path[0, :-1])
ax.set_xlabel("$t$")
ax.set_ylabel("$k_t$")
ax.set_title("Optimal path of $k_t$")
plt.show()

_images/4b43323bd819f1842fed2a3c8c9ae26ca510e0dd64829abfa27276cce3522202.png

Now we’ll depart from the preceding transition matrix that made the Markov state be strictly periodic.

We’ll begin with symmetric transition matrices of the form

\[\begin{split} \Pi_{2}=\begin{bmatrix} 1-\lambda & \lambda\\ \lambda & 1-\lambda \end{bmatrix}. \end{split}\]

λ = 0.8 # high λ
Π2 = np.array([[1-λ, λ],
               [λ, 1-λ]])

ex1_b = qe.LQMarkov(Π2, Qs, Rs, As, Bs, Cs=Cs, Ns=Ns, beta=β)
ex1_b.stationary_values();
ex1_b.Fs

array([[[ 0.57291724, -0.28645862]],

       [[ 0.74434525, -0.37217263]]])

λ = 0.2 # low λ
Π2 = np.array([[1-λ, λ],
               [λ, 1-λ]])

ex1_b = qe.LQMarkov(Π2, Qs, Rs, As, Bs, Cs=Cs, Ns=Ns, beta=β)
ex1_b.stationary_values();
ex1_b.Fs

array([[[ 0.59533259, -0.2976663 ]],

       [[ 0.72818728, -0.36409364]]])

We can plot optimal decision rules associated with different \(\lambda\) values.

λ_vals = np.linspace(0., 1., 10)
F1 = np.empty((λ_vals.size, 2))
F2 = np.empty((λ_vals.size, 2))

for i, λ in enumerate(λ_vals):
    Π2 = np.array([[1-λ, λ],
                   [λ, 1-λ]])

    ex1_b = qe.LQMarkov(Π2, Qs, Rs, As, Bs, Cs=Cs, Ns=Ns, beta=β)
    ex1_b.stationary_values();
    F1[i, :] = ex1_b.Fs[0, 0, :]
    F2[i, :] = ex1_b.Fs[1, 0, :]

for i, state_var in enumerate(state_vec1):
    fig, ax = plt.subplots()
    ax.plot(λ_vals, F1[:, i], label="$\overline{s}_1$", color="b")
    ax.plot(λ_vals, F2[:, i], label="$\overline{s}_2$", color="r")

    ax.set_xlabel("$\lambda$")
    ax.set_ylabel("$F_{s_t}$")
    ax.set_title(f"Coefficient on {state_var}")
    ax.legend()
    plt.show()

<>:3: SyntaxWarning: invalid escape sequence '\o'
<>:4: SyntaxWarning: invalid escape sequence '\o'
<>:6: SyntaxWarning: invalid escape sequence '\l'
<>:3: SyntaxWarning: invalid escape sequence '\o'
<>:4: SyntaxWarning: invalid escape sequence '\o'
<>:6: SyntaxWarning: invalid escape sequence '\l'
/tmp/ipykernel_7172/1661029301.py:3: SyntaxWarning: invalid escape sequence '\o'
  ax.plot(λ_vals, F1[:, i], label="$\overline{s}_1$", color="b")
/tmp/ipykernel_7172/1661029301.py:4: SyntaxWarning: invalid escape sequence '\o'
  ax.plot(λ_vals, F2[:, i], label="$\overline{s}_2$", color="r")
/tmp/ipykernel_7172/1661029301.py:6: SyntaxWarning: invalid escape sequence '\l'
  ax.set_xlabel("$\lambda$")

_images/134d58e5bb9a5d4e1a501f0b0e9104df807eb9232bf73275f91ee53a199a0da3.png

_images/17fe63c322edb97e4678b7500109320bc7d33c667bd3e9f546f3e531925dd653.png

Notice how the decision rules’ constants and slopes behave as functions of \(\lambda\).

Evidently, as the Markov chain becomes more nearly periodic (i.e., as \(\lambda \rightarrow 1\)), the dynamic program adjusts capital faster in the low adjustment cost Markov state to take advantage of what is only temporarily a more favorable time to invest.

Now let’s study situations in which the Markov transition matrix \(\Pi\) is asymmetric

\[\begin{split} \Pi_{3}=\begin{bmatrix} 1-\lambda & \lambda\\ \delta & 1-\delta \end{bmatrix}. \end{split}\]

λ, δ = 0.8, 0.2
Π3 = np.array([[1-λ, λ],
               [δ, 1-δ]])

ex1_b = qe.LQMarkov(Π3, Qs, Rs, As, Bs, Cs=Cs, Ns=Ns, beta=β)
ex1_b.stationary_values();
ex1_b.Fs

array([[[ 0.57169781, -0.2858489 ]],

       [[ 0.72749075, -0.36374537]]])

We can plot optimal decision rules for different \(\lambda\) and \(\delta\) values.

λ_vals = np.linspace(0., 1., 10)
δ_vals = np.linspace(0., 1., 10)

λ_grid = np.empty((λ_vals.size, δ_vals.size))
δ_grid = np.empty((λ_vals.size, δ_vals.size))
F1_grid = np.empty((λ_vals.size, δ_vals.size, len(state_vec1)))
F2_grid = np.empty((λ_vals.size, δ_vals.size, len(state_vec1)))

for i, λ in enumerate(λ_vals):
    λ_grid[i, :] = λ
    δ_grid[i, :] = δ_vals
    for j, δ in enumerate(δ_vals):
        Π3 = np.array([[1-λ, λ],
                       [δ, 1-δ]])

        ex1_b = qe.LQMarkov(Π3, Qs, Rs, As, Bs, Cs=Cs, Ns=Ns, beta=β)
        ex1_b.stationary_values();
        F1_grid[i, j, :] = ex1_b.Fs[0, 0, :]
        F2_grid[i, j, :] = ex1_b.Fs[1, 0, :]

for i, state_var in enumerate(state_vec1):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    # high adjustment cost, blue surface
    ax.plot_surface(λ_grid, δ_grid, F1_grid[:, :, i], color="b")
    # low adjustment cost, red surface
    ax.plot_surface(λ_grid, δ_grid, F2_grid[:, :, i], color="r")
    ax.set_xlabel("$\lambda$")
    ax.set_ylabel("$\delta$")
    ax.set_zlabel("$F_{s_t}$")
    ax.set_title(f"coefficient on {state_var}")
    plt.show()

<>:8: SyntaxWarning: invalid escape sequence '\l'
<>:9: SyntaxWarning: invalid escape sequence '\d'
<>:8: SyntaxWarning: invalid escape sequence '\l'
<>:9: SyntaxWarning: invalid escape sequence '\d'
/tmp/ipykernel_7172/3823440905.py:8: SyntaxWarning: invalid escape sequence '\l'
  ax.set_xlabel("$\lambda$")
/tmp/ipykernel_7172/3823440905.py:9: SyntaxWarning: invalid escape sequence '\d'
  ax.set_ylabel("$\delta$")

_images/6283cc4a9ca060dfd811bf113d5ac3ab8f42321239197c842fd9f5ac2eaa9578.png

_images/f62a64f689b13d1772f094d1f87edb4d4be70bad77bc89a933427bb267a1e4bc.png

The following code defines a wrapper function that computes optimal decision rules for cases with different Markov transition matrices

def run(construct_func, vals_dict, state_vec):
    """
    A Wrapper function that repeats the computation above
    for different cases
    """

    Qs, Rs, Ns, As, Bs, Cs, k_star = construct_func(**vals_dict)

    # Symmetric Π
    # Notice that pure periodic transition is a special case
    # when λ=1
    print("symmetric Π case:\n")
    λ_vals = np.linspace(0., 1., 10)
    F1 = np.empty((λ_vals.size, len(state_vec)))
    F2 = np.empty((λ_vals.size, len(state_vec)))

    for i, λ in enumerate(λ_vals):
        Π2 = np.array([[1-λ, λ],
                       [λ, 1-λ]])

        mplq = qe.LQMarkov(Π2, Qs, Rs, As, Bs, Cs=Cs, Ns=Ns, beta=β)
        mplq.stationary_values();
        F1[i, :] = mplq.Fs[0, 0, :]
        F2[i, :] = mplq.Fs[1, 0, :]

    for i, state_var in enumerate(state_vec):
        fig = plt.figure()
        ax = fig.add_subplot(111)
        ax.plot(λ_vals, F1[:, i], label="$\overline{s}_1$", color="b")
        ax.plot(λ_vals, F2[:, i], label="$\overline{s}_2$", color="r")

        ax.set_xlabel("$\lambda$")
        ax.set_ylabel("$F(\overline{s}_t)$")
        ax.set_title(f"coefficient on {state_var}")
        ax.legend()
        plt.show()

    # Plot optimal k*_{s_t} and k that optimal policies are targeting
    # only for example 1
    if state_vec == ["k", "constant term"]:
        fig = plt.figure()
        ax = fig.add_subplot(111)
        for i in range(2):
            F = [F1, F2][i]
            c = ["b", "r"][i]
            ax.plot([0, 1], [k_star[i], k_star[i]], "--",
                    color=c, label="$k^*(\overline{s}_"+str(i+1)+")$")
            ax.plot(λ_vals, - F[:, 1] / F[:, 0], color=c,
                    label="$k^{target}(\overline{s}_"+str(i+1)+")$")

        # Plot a vertical line at λ=0.5
        ax.plot([0.5, 0.5], [min(k_star), max(k_star)], "-.")

        ax.set_xlabel("$\lambda$")
        ax.set_ylabel("$k$")
        ax.set_title("Optimal k levels and k targets")
        ax.text(0.5, min(k_star)+(max(k_star)-min(k_star))/20, "$\lambda=0.5$")
        ax.legend(bbox_to_anchor=(1., 1.))
        plt.show()

    # Asymmetric Π
    print("asymmetric Π case:\n")
    δ_vals = np.linspace(0., 1., 10)

    λ_grid = np.empty((λ_vals.size, δ_vals.size))
    δ_grid = np.empty((λ_vals.size, δ_vals.size))
    F1_grid = np.empty((λ_vals.size, δ_vals.size, len(state_vec)))
    F2_grid = np.empty((λ_vals.size, δ_vals.size, len(state_vec)))

    for i, λ in enumerate(λ_vals):
        λ_grid[i, :] = λ
        δ_grid[i, :] = δ_vals
        for j, δ in enumerate(δ_vals):
            Π3 = np.array([[1-λ, λ],
                           [δ, 1-δ]])

            mplq = qe.LQMarkov(Π3, Qs, Rs, As, Bs, Cs=Cs, Ns=Ns, beta=β)
            mplq.stationary_values();
            F1_grid[i, j, :] = mplq.Fs[0, 0, :]
            F2_grid[i, j, :] = mplq.Fs[1, 0, :]

    for i, state_var in enumerate(state_vec):
        fig = plt.figure()
        ax = fig.add_subplot(111, projection='3d')
        ax.plot_surface(λ_grid, δ_grid, F1_grid[:, :, i], color="b")
        ax.plot_surface(λ_grid, δ_grid, F2_grid[:, :, i], color="r")
        ax.set_xlabel("$\lambda$")
        ax.set_ylabel("$\delta$")
        ax.set_zlabel("$F(\overline{s}_t)$")
        ax.set_title(f"coefficient on {state_var}")
        plt.show()

<>:29: SyntaxWarning: invalid escape sequence '\o'
<>:30: SyntaxWarning: invalid escape sequence '\o'
<>:32: SyntaxWarning: invalid escape sequence '\l'
<>:33: SyntaxWarning: invalid escape sequence '\o'
<>:47: SyntaxWarning: invalid escape sequence '\o'
<>:49: SyntaxWarning: invalid escape sequence '\o'
<>:54: SyntaxWarning: invalid escape sequence '\l'
<>:57: SyntaxWarning: invalid escape sequence '\l'
<>:87: SyntaxWarning: invalid escape sequence '\l'
<>:88: SyntaxWarning: invalid escape sequence '\d'
<>:89: SyntaxWarning: invalid escape sequence '\o'
<>:29: SyntaxWarning: invalid escape sequence '\o'
<>:30: SyntaxWarning: invalid escape sequence '\o'
<>:32: SyntaxWarning: invalid escape sequence '\l'
<>:33: SyntaxWarning: invalid escape sequence '\o'
<>:47: SyntaxWarning: invalid escape sequence '\o'
<>:49: SyntaxWarning: invalid escape sequence '\o'
<>:54: SyntaxWarning: invalid escape sequence '\l'
<>:57: SyntaxWarning: invalid escape sequence '\l'
<>:87: SyntaxWarning: invalid escape sequence '\l'
<>:88: SyntaxWarning: invalid escape sequence '\d'
<>:89: SyntaxWarning: invalid escape sequence '\o'
/tmp/ipykernel_7172/3965109342.py:29: SyntaxWarning: invalid escape sequence '\o'
  ax.plot(λ_vals, F1[:, i], label="$\overline{s}_1$", color="b")
/tmp/ipykernel_7172/3965109342.py:30: SyntaxWarning: invalid escape sequence '\o'
  ax.plot(λ_vals, F2[:, i], label="$\overline{s}_2$", color="r")
/tmp/ipykernel_7172/3965109342.py:32: SyntaxWarning: invalid escape sequence '\l'
  ax.set_xlabel("$\lambda$")
/tmp/ipykernel_7172/3965109342.py:33: SyntaxWarning: invalid escape sequence '\o'
  ax.set_ylabel("$F(\overline{s}_t)$")
/tmp/ipykernel_7172/3965109342.py:47: SyntaxWarning: invalid escape sequence '\o'
  color=c, label="$k^*(\overline{s}_"+str(i+1)+")$")
/tmp/ipykernel_7172/3965109342.py:49: SyntaxWarning: invalid escape sequence '\o'
  label="$k^{target}(\overline{s}_"+str(i+1)+")$")
/tmp/ipykernel_7172/3965109342.py:54: SyntaxWarning: invalid escape sequence '\l'
  ax.set_xlabel("$\lambda$")
/tmp/ipykernel_7172/3965109342.py:57: SyntaxWarning: invalid escape sequence '\l'
  ax.text(0.5, min(k_star)+(max(k_star)-min(k_star))/20, "$\lambda=0.5$")
/tmp/ipykernel_7172/3965109342.py:87: SyntaxWarning: invalid escape sequence '\l'
  ax.set_xlabel("$\lambda$")
/tmp/ipykernel_7172/3965109342.py:88: SyntaxWarning: invalid escape sequence '\d'
  ax.set_ylabel("$\delta$")
/tmp/ipykernel_7172/3965109342.py:89: SyntaxWarning: invalid escape sequence '\o'
  ax.set_zlabel("$F(\overline{s}_t)$")

To illustrate the code with another example, we shall set \(f_{2,s_t}\) and \(d_{s_t}\) as constant functions and

\[ f_{1,1} = 0.5, f_{1,2} = 1 \]

Thus, the sole role of the Markov jump state \(s_t\) is to identify times in which capital is very productive and other times in which it is less productive.

The example below reveals much about the structure of the optimum problem and optimal policies.

Only \(f_{1,{s_t}}\) varies with \(s_t\).

So there are different \(s_t\)-dependent optimal static \(k\) level in different states \(k^*_{s_t}=\frac{f_{1,{s_t}}}{2 f_{2,{s_t}}}\), values of \(k\) that maximize one-period payoff functions in each state.

We denote a target \(k\) level as \(k^{target}_{s_t}\), the fixed point of the optimal policies in each state, given the value of \(\lambda\).

We call \(k^{target}_{s_t}\) a “target” because in each Markov state \(s_t\), optimal policies are contraction mappings and will push \(k_t\) towards a fixed point \(k^{target}_{s_t}\).

When \(\lambda \rightarrow 0\), each Markov state becomes close to absorbing state and consequently \(k^{target}_{s_t} \rightarrow k^*_{s_t}\).

But when \(\lambda \rightarrow 1\), the Markov transition matrix becomes more nearly periodic, so the optimum decision rules target more at the optimal \(k\) level in the other state in order to enjoy higher expected payoff in the next period.

The switch happens at \(\lambda = 0.5\) when both states are equally likely to be reached.

Below we plot an additional figure that shows optimal \(k\) levels in the two states Markov jump state and also how the targeted \(k\) levels change as \(\lambda\) changes.

run(construct_arrays1, {"f1_vals":[0.5, 1.]}, state_vec1)

symmetric Π case:

_images/b27302ccbfbd9b3c2c30dca0197c9ef373e4132914480c1628dbeab9070fc649.png

_images/4dc9976641fdc3614660167aee20db748d554ada34e314bdd6d3a85645db48b6.png

_images/d1446a9d7118e84152a00f3cd3666555ee2453229063ebd9ae9805933ae87a7a.png

asymmetric Π case:

_images/a6c194e104f2199ea87c336a5830633b2052f466fd291411c255aa27f0df3064.png

_images/aad447205db80b20a3c91201fc96f405b17972cfa0add1dbdd842120ac2648d3.png

Set \(f_{1,{s_t}}\) and \(d_{s_t}\) as constant functions and

\[ f_{2,1} = 0.5, f_{2,2} = 1 \]

run(construct_arrays1, {"f2_vals":[0.5, 1.]}, state_vec1)

symmetric Π case:

_images/dba8506d886eebdbaca2804f9a27d3d282bf21f1f2edcdbf4d0ce56caeb59f65.png

_images/87d9218fbc30cf597d81fd8669c63d9ce504a521bf18810dca6472211fff683d.png

_images/ce43701af54f120543485a5df76b27e5fef60713ce711e28279031856f949083.png

asymmetric Π case:

_images/97618bf0b5d4fdc60d634572ae50ce7946e664bdb2da50245e01ec13c82f534e.png

_images/f3ba1df1aba60178b316ff5a3c26c953c85a1a33c3f7120a12ea0d0e0e8ff0c8.png

8.6. Example 2#

We now add to the example 1 setup another state variable \(w_t\) that follows the evolution law

\[ w_{t+1}=\alpha_{0}\left(s_{t}\right)+\rho\left(s_{t}\right)w_{t}+\sigma\left(s_{t}\right)\epsilon_{t+1},\quad\epsilon_{t+1}\sim N\left(0,1\right) \]

We think of \(w_t\) as a rental rate or tax rate that the decision maker pays each period for \(k_t\).

To capture this idea, we add to the decision-maker’s one-period payoff function the product of \(w_t\) and \(k_t\)

\[ r(s_{t},k_{t},w_{t})=f_{1,s_{t}} k_{t}-f_{2, s_{t}} k_{t}^{2}-d_{s_{t}} (k_{t+1}-k_{t})^{2}-w_{t}k_{t}, \]

We now let the continuous part of the state at time \(t\) be \(x_{t}=\begin{bmatrix} k_{t}\\ 1\\ w_{t} \end{bmatrix}\) and continue to set the control \(u_{t}=k_{t+1}-k_{t}\).

We can write the one-period payoff function \(r\left(s_{t},k_{t},w_{t}\right)\) and the state-transition law as

\[\begin{split} \begin{aligned} r\left(s_{t},k_{t},w_{t}\right) & =f_{1}\left(s_{t}\right)k_{t}-f_{2}\left(s_{t}\right)k_{t}^{2}-d\left(s_{t}\right)\left(k_{t+1}-k_{t}\right)^{2}-w_{t}k_{t} \\ & =-\left(x_{t}^{\prime}\underset{\equiv R\left(s_{t}\right)}{\underbrace{ \begin{bmatrix} f_{2}\left(s_{t}\right) & -\frac{f_{1}\left(s_{t}\right)}{2} & \frac{1}{2}\\ -\frac{f_{1}\left(s_{t}\right)}{2} & 0 & 0\\ \frac{1}{2} & 0 & 0 \end{bmatrix}}} x_{t}+ \underset{\equiv Q\left(s_{t}\right)}{\underbrace{d\left(s_{t}\right)}}u_{t}^{2}\right), \end{aligned} \end{split}\]

and

\[\begin{split} x_{t+1}=\begin{bmatrix} k_{t+1}\\ 1\\ w_{t+1} \end{bmatrix}=\underset{\equiv A\left(s_{t}\right)}{\underbrace{\begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & \alpha_{0}\left(s_{t}\right) & \rho\left(s_{t}\right) \end{bmatrix}}}x_{t}+\underset{\equiv B\left(s_{t}\right)}{\underbrace{\begin{bmatrix} 1\\ 0\\ 0 \end{bmatrix}}}u_{t}+\underset{\equiv C\left(s_{t}\right)}{\underbrace{\begin{bmatrix} 0\\ 0\\ \sigma\left(s_{t}\right) \end{bmatrix}}}\epsilon_{t+1} \end{split}\]

def construct_arrays2(f1_vals=[1. ,1.],
                      f2_vals=[1., 1.],
                      d_vals=[1., 1.],
                      α0_vals=[1., 1.],
                      ρ_vals=[0.9, 0.9],
                      σ_vals=[1., 1.]):
    """
    Construct matrices that maps the problem described in example 2
    into a Markov jump linear quadratic dynamic programming problem.
    """

    m = len(f1_vals)
    n, k, j = 3, 1, 1

    Rs = np.zeros((m, n, n))
    Qs = np.zeros((m, k, k))
    As = np.zeros((m, n, n))
    Bs = np.zeros((m, n, k))
    Cs = np.zeros((m, n, j))

    for i in range(m):
        Rs[i, 0, 0] = f2_vals[i]
        Rs[i, 1, 0] = - f1_vals[i] / 2
        Rs[i, 0, 1] = - f1_vals[i] / 2
        Rs[i, 0, 2] = 1/2
        Rs[i, 2, 0] = 1/2

        Qs[i, 0, 0] = d_vals[i]

        As[i, 0, 0] = 1
        As[i, 1, 1] = 1
        As[i, 2, 1] = α0_vals[i]
        As[i, 2, 2] = ρ_vals[i]

        Bs[i, :, :] = np.array([[1, 0, 0]]).T

        Cs[i, :, :] = np.array([[0, 0, σ_vals[i]]]).T

    Ns = None
    k_star = None

    return Qs, Rs, Ns, As, Bs, Cs, k_star

state_vec2 = ["k", "constant term", "w"]

Only \(d_{s_t}\) depends on \(s_t\).

run(construct_arrays2, {"d_vals":[1., 0.5]}, state_vec2)

symmetric Π case:

_images/b38c41a68b4f2b7efaefa351e65d62b28c6b64eafdeb07a75c030003951e6255.png

_images/263fd36b53910623f3763145aabea5b09e3d78e22f25b208762b7ed1baab01a8.png

_images/34fc0d00e00900c41966112a73bc68c4fb89946c5fd8c9b4673f6df31296da04.png

asymmetric Π case:

_images/fd2be1a2249a3baf1e1925a2618eb8a3f0b18a18a685a61f132d06269991d033.png

_images/e2e74516253d6245967ebf28899238cbb7882e7a3b005a83eb6ea67d248d6fea.png

_images/be8fa9b038e287a8e46c7e3078da5dae0558ad0cb6b7ea32ea8a8fa8672d67b2.png

Only \(f_{1,{s_t}}\) depends on \(s_t\).

run(construct_arrays2, {"f1_vals":[0.5, 1.]}, state_vec2)

symmetric Π case:

_images/6623b698ed7c60f6a7e833061a2a8c574f22d44fca5cc88e3f7f5d18d339d870.png

_images/436330f23fee2249430c94612fcc4a0959e66c6fa0d35b8f58a8d836fdbbc9b4.png

asymmetric Π case:

_images/96a08070383d7c23d8e6a5adebe053ea540e3eecc016a8909133e36b14d350d0.png

_images/70576603fb07f750502e72a867e8f22e6c71dde816f33d89a3e0fcbdf1642329.png

Only \(f_{2,{s_t}}\) depends on \(s_t\).

run(construct_arrays2, {"f2_vals":[0.5, 1.]}, state_vec2)

symmetric Π case:

_images/61b6b3243edb2f66e4b1b1db73ada2dc650b2ab36ac298e89cd53ef3aef21f27.png

_images/418f2f210ea81a63ee2d6d08f8b2137808060141177cc197cb3fb223f690e561.png

asymmetric Π case:

_images/28413c4e3cf2f88e5442d773720bef5b7f884d4c8dc7fee4d5a56eb72133a23e.png

_images/008932b152e723a7bf41575f5f980f7ef4f70c9475211237a170ed6229fd9738.png

Only \(\alpha_0(s_t)\) depends on \(s_t\).

run(construct_arrays2, {"α0_vals":[0.5, 1.]}, state_vec2)

symmetric Π case:

_images/cb118b38b366ff4453aa0a5aabca56c8d7fad1dac3dbec307758ad9dbed810b6.png

asymmetric Π case:

_images/22c8e6ae0b4728127582e6f9e67b1974031cb057ed7a0c88c3d0a0559c3e9d0c.png

Only \(\rho_{s_t}\) depends on \(s_t\).

run(construct_arrays2, {"ρ_vals":[0.5, 0.9]}, state_vec2)

symmetric Π case:

_images/66c1951774d2b4a806cb21091b982e66319c919772d366cf2cb6e2de136eba85.png

_images/6e69423146ebd7990c38a41c6fa048a1a1613ef7f8ce0947bceea83ef0a6750d.png

asymmetric Π case:

_images/61031f2ca7e4e08d9745c6fe5f439536f9de7319033fb024b12b65901f84370b.png

_images/ce8d8f11aa089ce7178f768b4a5cfd0bcc6e3d8bde0d61f65645085ab4c061f1.png

Only \(\sigma_{s_t}\) depends on \(s_t\).

run(construct_arrays2, {"σ_vals":[0.5, 1.]}, state_vec2)

symmetric Π case:

_images/e684721b6068f43b2a5afa2f8065460074662cb87ac17b70e2a760bb379e02c7.png

asymmetric Π case:

_images/78b4a4c28ccab10679a0a734a20cdbdc3d3a9c84663120b8ba5ddc6b137dc12e.png

8.7. More examples#

The following lectures describe how Markov jump linear quadratic dynamic programming can be used to extend the [Barro, 1979] model of optimal tax-smoothing and government debt in several interesting directions