diff --git a/lectures/lqcontrol.md b/lectures/lqcontrol.md index 79c2f800e..f04dff5c4 100644 --- a/lectures/lqcontrol.md +++ b/lectures/lqcontrol.md @@ -50,14 +50,14 @@ These themes appear repeatedly below. Mathematically, LQ control problems are closely related to {doc}`the Kalman filter ` -* Recursive formulations of linear-quadratic control problems and Kalman filtering problems both involve matrix [Riccati equations](https://en.wikipedia.org/wiki/Riccati_equation). -* Classical formulations of linear control and linear filtering problems make use of similar matrix decompositions (see for example[Classical Control with Linear Algebra](https://python-advanced.quantecon.org/lu_tricks.html) and [Classical Prediction and Filtering With Linear Algebra](https://python-advanced.quantecon.org/classical_filtering.html)). +* Recursive formulations of linear-quadratic control problems and Kalman filtering problems both involve matrix **Riccati equations**. +* Classical formulations of linear control and linear filtering problems make use of similar matrix decompositions (see for example [this lecture](https://python-advanced.quantecon.org/lu_tricks.html) and [this lecture](https://python-advanced.quantecon.org/classical_filtering.html)). In reading what follows, it will be useful to have some familiarity with * matrix manipulations * vectors of random variables -* dynamic programming and the Bellman equation (see for example {doc}`Shortest Paths ` and {doc}`Optimal Growth `) +* dynamic programming and the Bellman equation (see for example {doc}`this lecture ` and {doc}`this lecture `) For additional reading on LQ control, see, for example, @@ -97,7 +97,7 @@ x_{t+1} = A x_t + B u_t + C w_{t+1}, Here * $u_t$ is a "control" vector, incorporating choices available to a decision-maker confronting the current state $x_t$ -* $\{w_t\}$ is a sequence of uncorrelated zero mean shock process satisfying $\mathbb E w_t w_t^\top = I$, where the right-hand side is the identity matrix +* $\{w_t\}$ is an uncorrelated zero mean shock process satisfying $\mathbb E w_t w_t' = I$, where the right-hand side is the identity matrix Regarding the dimensions @@ -116,8 +116,8 @@ $$ Here $a_t$ is assets, $r$ is a fixed interest rate, $c_t$ is current consumption, and $y_t$ is current non-financial income. -If we suppose that $\{y_t\}$ is serially uncorrelated and $N(0, -\sigma^2)$, then, taking $\{w_t\}$ to be standard normal, we can write +If we suppose that $\{ y_t \}$ is serially uncorrelated and $N(0, +\sigma^2)$, then, taking $\{ w_t \}$ to be standard normal, we can write the system as $$ @@ -250,11 +250,9 @@ In the LQ model, the aim is to minimize flow of losses, where time-$t$ loss is g ```{math} :label: lq_pref_flow -x_t^\top R x_t + u_t^\top Q u_t +x_t' R x_t + u_t' Q u_t ``` -Where the entries in matrices $R$ and $Q$ are chosen to match the specific problem that is being studied. - Here * $R$ is assumed to be $n \times n$, symmetric and nonnegative definite. @@ -262,11 +260,6 @@ Here ```{note} In fact, for many economic problems, the definiteness conditions on $R$ and $Q$ can be relaxed. It is sufficient that certain submatrices of $R$ and $Q$ be nonnegative definite. See {cite}`HansenSargent2008` for details. - -``` - -```{note} -The use of $R$ and $Q$ notations may differ depending on the source. Some authors use $Q$ to be the matrix associated with the state variables and $R$ with the control variables. ``` #### Example 1 @@ -275,7 +268,7 @@ A very simple example that satisfies these assumptions is to take $R$ and $Q$ to be identity matrices so that current loss is $$ -x_t^\top I x_t + u_t^\top I u_t = \| x_t \|^2 + \| u_t \|^2 +x_t' I x_t + u_t' I u_t = \| x_t \|^2 + \| u_t \|^2 $$ Thus, for both the state and the control, loss is measured as squared distance from the origin. @@ -297,7 +290,7 @@ In the household problem {ref}`studied above `, setting $R=0$ and $Q=1$ yields preferences $$ -x_t^\top R x_t + u_t^\top Q u_t = u_t^2 = (c_t - \bar c)^2 +x_t' R x_t + u_t' Q u_t = u_t^2 = (c_t - \bar c)^2 $$ Under this specification, the household's current loss is the squared deviation of consumption from the ideal level $\bar c$. @@ -320,7 +313,7 @@ In this case, the aim is to choose a sequence of controls $\{u_0, \ldots, u_{T-1 \mathbb E \, \left\{ - \sum_{t=0}^{T-1} \beta^t (x_t^\top R x_t + u_t^\top Q u_t) + \beta^T x_T^\top R_f x_T + \sum_{t=0}^{T-1} \beta^t (x_t' R x_t + u_t' Q u_t) + \beta^T x_T' R_f x_T \right\} ``` @@ -328,13 +321,13 @@ subject to the law of motion {eq}`lq_lom` and initial state $x_0$. The new objects introduced here are $\beta$ and the matrix $R_f$. -The scalar $\beta$ is the discount factor, while $x^\top R_f x$ gives terminal loss associated with state $x$. +The scalar $\beta$ is the discount factor, while $x' R_f x$ gives terminal loss associated with state $x$. Comments: * We assume $R_f$ to be $n \times n$, symmetric and nonnegative definite. * We allow $\beta = 1$, and hence include the undiscounted case. -* $x_0$ may itself be random, in which case we require it to be independent of the shock sequence $\{w_1, \ldots, w_T\}$. +* $x_0$ may itself be random, in which case we require it to be independent of the shock sequence $w_1, \ldots, w_T$. (lq_cp)= ### Information @@ -375,9 +368,9 @@ What's special about the LQ case is that -- as we shall soon see --- the optima ### Solution To solve the finite horizon LQ problem we can use a dynamic programming -strategy based on backward induction that is conceptually similar to the approach adopted in {doc}`Shortest Paths `. +strategy based on backward induction that is conceptually similar to the approach adopted in {doc}`this lecture `. -For reasons that will soon become clear, we first introduce the notation $J_T(x) = x^\top R_f x$. +For reasons that will soon become clear, we first introduce the notation $J_T(x) = x' R_f x$. Now consider the problem of the decision-maker in the second to last period. @@ -389,13 +382,11 @@ solves $$ \min_u \{ -x_{T-1}^\top R x_{T-1} + u^\top Q u + \beta \, +x_{T-1}' R x_{T-1} + u' Q u + \beta \, \mathbb E J_T(A x_{T-1} + B u + C w_T) \} $$ -We use $u$ instead of $u_{T-1}$ here to simplify notation. - At this stage, it is convenient to define the function ```{math} @@ -403,7 +394,7 @@ At this stage, it is convenient to define the function J_{T-1} (x) = \min_u \{ -x^\top R x + u^\top Q u + \beta \, +x' R x + u' Q u + \beta \, \mathbb E J_T(A x + B u + C w_T) \} ``` @@ -412,7 +403,7 @@ The function $J_{T-1}$ will be called the $T-1$ value function, and $J_{T-1}(x)$ Now let's step back to $T-2$. -For a decision-maker at $T-2$, the value $J_{T-1}(x)$ plays a role analogous to that played by the terminal loss $J_T(x) = x^\top R_f x$ for the decision-maker at $T-1$. +For a decision-maker at $T-2$, the value $J_{T-1}(x)$ plays a role analogous to that played by the terminal loss $J_T(x) = x' R_f x$ for the decision-maker at $T-1$. That is, $J_{T-1}(x)$ summarizes the future loss associated with moving to state $x$. @@ -426,7 +417,7 @@ Her problem is therefore $$ \min_u \{ -x_{T-2}^\top R x_{T-2} + u^\top Q u + \beta \, +x_{T-2}' R x_{T-2} + u' Q u + \beta \, \mathbb E J_{T-1}(Ax_{T-2} + B u + C w_{T-1}) \} $$ @@ -437,7 +428,7 @@ $$ J_{T-2} (x) = \min_u \{ -x^\top R x + u^\top Q u + \beta \, +x' R x + u' Q u + \beta \, \mathbb E J_{T-1}(Ax + B u + C w_{T-1}) \} $$ @@ -450,11 +441,11 @@ $$ J_{t-1} (x) = \min_u \{ -x^\top R x + u^\top Q u + \beta \, +x' R x + u' Q u + \beta \, \mathbb E J_{t}(Ax + B u + C w_t) \} \quad \text{and} \quad -J_T(x) = x^\top R_f x +J_T(x) = x' R_f x $$ The first equality is the Bellman equation from dynamic programming theory specialized to the finite horizon LQ problem. @@ -463,7 +454,7 @@ Now that we have $\{J_0, \ldots, J_T\}$, we can obtain the optimal controls. As a first step, let's find out what the value functions look like. -It turns out that every $J_t$ has the form $J_t(x) = x^\top P_t x + d_t$ where $P_t$ is a $n \times n$ matrix and $d_t$ is a constant. +It turns out that every $J_t$ has the form $J_t(x) = x' P_t x + d_t$ where $P_t$ is a $n \times n$ matrix and $d_t$ is a constant. We can show this by induction, starting from $P_T := R_f$ and $d_T = 0$. @@ -474,8 +465,8 @@ Using this notation, {eq}`lq_lsm` becomes J_{T-1} (x) = \min_u \{ -x^\top R x + u^\top Q u + \beta \, -\mathbb E (A x + B u + C w_T)^\top P_T (A x + B u + C w_T) +x' R x + u' Q u + \beta \, +\mathbb E (A x + B u + C w_T)' P_T (A x + B u + C w_T) \} ``` @@ -486,13 +477,13 @@ Applying the relevant rules of {ref}`matrix calculus `, this gives ```{math} :label: lq_oc0 -u = - (Q + \beta B^\top P_T B)^{-1} \beta B^\top P_T A x +u = - (Q + \beta B' P_T B)^{-1} \beta B' P_T A x ``` Plugging this back into {eq}`lq_fswb` and rearranging yields $$ -J_{T-1} (x) = x^\top P_{T-1} x + d_{T-1} +J_{T-1} (x) = x' P_{T-1} x + d_{T-1} $$ where @@ -500,8 +491,8 @@ where ```{math} :label: lq_finr -P_{T-1} = R - \beta^2 A^\top P_T B (Q + \beta B^\top P_T B)^{-1} B^\top P_T A + -\beta A^\top P_T A +P_{T-1} = R - \beta^2 A' P_T B (Q + \beta B' P_T B)^{-1} B' P_T A + +\beta A' P_T A ``` and @@ -509,18 +500,18 @@ and ```{math} :label: lq_finrd -d_{T-1} := \beta \mathop{\mathrm{trace}}(C^\top P_T C) +d_{T-1} := \beta \mathop{\mathrm{trace}}(C' P_T C) ``` (The algebra is a good exercise --- we'll leave it up to you.) -If we continue working backwards in this manner, it soon becomes clear that $J_t (x) = x^\top P_t x + d_t$ as claimed, where $P_t$ and $d_t$ satisfy the recursions +If we continue working backwards in this manner, it soon becomes clear that $J_t (x) = x' P_t x + d_t$ as claimed, where $\{P_t\}$ and $\{d_t\}$ satisfy the recursions ```{math} :label: lq_pr -P_{t-1} = R - \beta^2 A^\top P_t B (Q + \beta B^\top P_t B)^{-1} B^\top P_t A + -\beta A^\top P_t A +P_{t-1} = R - \beta^2 A' P_t B (Q + \beta B' P_t B)^{-1} B' P_t A + +\beta A' P_t A \quad \text{with } \quad P_T = R_f ``` @@ -530,7 +521,7 @@ and ```{math} :label: lq_dd -d_{t-1} = \beta (d_t + \mathop{\mathrm{trace}}(C^\top P_t C)) +d_{t-1} = \beta (d_t + \mathop{\mathrm{trace}}(C' P_t C)) \quad \text{with } \quad d_T = 0 ``` @@ -542,14 +533,14 @@ Recalling {eq}`lq_oc0`, the minimizers from these backward steps are u_t = - F_t x_t \quad \text{where} \quad -F_t := (Q + \beta B^\top P_{t+1} B)^{-1} \beta B^\top P_{t+1} A +F_t := (Q + \beta B' P_{t+1} B)^{-1} \beta B' P_{t+1} A ``` These are the linear optimal control policies we {ref}`discussed above `. In particular, the sequence of controls given by {eq}`lq_oc` and {eq}`lq_lom` solves our finite horizon LQ problem. -Rephrasing this more precisely, the sequence $\{u_0, \ldots, u_{T-1}\}$ given by +Rephrasing this more precisely, the sequence $u_0, \ldots, u_{T-1}$ given by ```{math} :label: lq_xud @@ -619,7 +610,7 @@ Here $q$ is a large positive constant, the role of which is to induce the consum As before we set $y_t = \sigma w_{t+1} + \mu$ and $u_t := c_t - \bar c$, after which the constraint can be written as in {eq}`lq_lomwc`. We saw how this constraint could be manipulated into the LQ formulation $x_{t+1} = -Ax_t + Bu_t + Cw_{t+1}$ by setting $x_t = (a_t \; 1)^\top$ and using the definitions in {eq}`lq_lowmc2`. +Ax_t + Bu_t + Cw_{t+1}$ by setting $x_t = (a_t \; 1)'$ and using the definitions in {eq}`lq_lowmc2`. To match with this state and control, the objective function {eq}`lq_pio` can be written in the form of {eq}`lq_object` by choosing @@ -628,22 +619,26 @@ $$ Q := 1, \quad R := -\begin{bmatrix} +\left( +\begin{array}{cc} 0 & 0 \\ 0 & 0 -\end{bmatrix}, +\end{array} +\right), \quad \text{and} \quad R_f := -\begin{bmatrix} +\left( +\begin{array}{cc} q & 0 \\ 0 & 0 -\end{bmatrix} +\end{array} +\right) $$ Now that the problem is expressed in LQ form, we can proceed to the solution by applying {eq}`lq_pr` and {eq}`lq_oc`. -After generating shocks $\{w_1, \ldots, w_T\}$, the dynamics for assets and +After generating shocks $w_1, \ldots, w_T$, the dynamics for assets and consumption can be simulated via {eq}`lq_xud`. The following figure was computed using $r = 0.05, \beta = 1 / (1+ r), @@ -654,7 +649,7 @@ The shocks $\{w_t\}$ were taken to be IID and standard normal. ```{code-cell} python3 # Model parameters r = 0.05 -β = 1 / (1 + r) +β = 1/(1 + r) T = 45 c_bar = 2 σ = 0.25 @@ -693,14 +688,14 @@ bbox = (0., 1.02, 1., .102) legend_args = {'bbox_to_anchor': bbox, 'loc': 3, 'mode': 'expand'} p_args = {'lw': 2, 'alpha': 0.7} -axes[0].plot(range(1, T+1), income, 'g-', label="non-financial income", +axes[0].plot(list(range(1, T+1)), income, 'g-', label="non-financial income", **p_args) -axes[0].plot(range(T), c, 'k-', label="consumption", **p_args) +axes[0].plot(list(range(T)), c, 'k-', label="consumption", **p_args) -axes[1].plot(range(1, T+1), np.cumsum(income - μ), 'r-', +axes[1].plot(list(range(1, T+1)), np.cumsum(income - μ), 'r-', label="cumulative unanticipated income", **p_args) -axes[1].plot(range(T+1), assets, 'b-', label="assets", **p_args) -axes[1].plot(range(T), np.zeros(T), 'k-') +axes[1].plot(list(range(T+1)), assets, 'b-', label="assets", **p_args) +axes[1].plot(list(range(T)), np.zeros(T), 'k-') for ax in axes: ax.grid() @@ -761,14 +756,14 @@ bbox = (0., 1.02, 1., .102) legend_args = {'bbox_to_anchor': bbox, 'loc': 3, 'mode': 'expand'} p_args = {'lw': 2, 'alpha': 0.7} -axes[0].plot(range(1, T+1), income, 'g-', label="non-financial income", +axes[0].plot(list(range(1, T+1)), income, 'g-', label="non-financial income", **p_args) -axes[0].plot(range(T), c, 'k-', label="consumption", **p_args) +axes[0].plot(list(range(T)), c, 'k-', label="consumption", **p_args) -axes[1].plot(range(1, T+1), np.cumsum(income - μ), 'r-', +axes[1].plot(list(range(1, T+1)), np.cumsum(income - μ), 'r-', label="cumulative unanticipated income", **p_args) -axes[1].plot(range(T+1), assets, 'b-', label="assets", **p_args) -axes[1].plot(range(T), np.zeros(T), 'k-') +axes[1].plot(list(range(T+1)), assets, 'b-', label="assets", **p_args) +axes[1].plot(list(range(T)), np.zeros(T), 'k-') for ax in axes: ax.grid() @@ -804,14 +799,14 @@ For further examples and a more systematic treatment, see {cite}`HansenSargent20 (lq_cpt)= ### Adding a Cross-Product Term -In some LQ problems, preferences include a cross-product term $u_t^\top N x_t$, so that the objective function becomes +In some LQ problems, preferences include a cross-product term $u_t' N x_t$, so that the objective function becomes ```{math} :label: lq_object_cp \mathbb E \, \left\{ - \sum_{t=0}^{T-1} \beta^t (x_t^\top R x_t + u_t^\top Q u_t + 2 u_t^\top N x_t) + \beta^T x_T^\top R_f x_T + \sum_{t=0}^{T-1} \beta^t (x_t' R x_t + u_t' Q u_t + 2 u_t' N x_t) + \beta^T x_T' R_f x_T \right\} ``` @@ -822,9 +817,9 @@ The sequence $\{P_t\}$ from {eq}`lq_pr` becomes ```{math} :label: lq_pr_cp -P_{t-1} = R - (\beta B^\top P_t A + N)^\top -(Q + \beta B^\top P_t B)^{-1} (\beta B^\top P_t A + N) + -\beta A^\top P_t A +P_{t-1} = R - (\beta B' P_t A + N)' +(Q + \beta B' P_t B)^{-1} (\beta B' P_t A + N) + +\beta A' P_t A \quad \text{with } \quad P_T = R_f ``` @@ -836,7 +831,7 @@ The policies in {eq}`lq_oc` are modified to u_t = - F_t x_t \quad \text{where} \quad -F_t := (Q + \beta B^\top P_{t+1} B)^{-1} (\beta B^\top P_{t+1} A + N) +F_t := (Q + \beta B' P_{t+1} B)^{-1} (\beta B' P_{t+1} A + N) ``` The sequence $\{d_t\}$ is unchanged from {eq}`lq_dd`. @@ -857,7 +852,7 @@ objective function given by \mathbb E \, \left\{ - \sum_{t=0}^{\infty} \beta^t (x_t^\top R x_t + u_t^\top Q u_t + 2 u_t^\top N x_t) + \sum_{t=0}^{\infty} \beta^t (x_t' R x_t + u_t' Q u_t + 2 u_t' N x_t) \right\} ``` @@ -879,9 +874,9 @@ The stationary matrix $P$ is the solution to the ```{math} :label: lq_pr_ih -P = R - (\beta B^\top P A + N)^\top -(Q + \beta B^\top P B)^{-1} (\beta B^\top P A + N) + -\beta A^\top P A +P = R - (\beta B' P A + N)' +(Q + \beta B' P B)^{-1} (\beta B' P A + N) + +\beta A' P A ``` Equation {eq}`lq_pr_ih` is also called the *LQ Bellman equation*, and the map @@ -895,7 +890,7 @@ The stationary optimal policy for this model is u = - F x \quad \text{where} \quad -F = (Q + \beta B^\top P B)^{-1} (\beta B^\top P A + N) +F = (Q + \beta B' P B)^{-1} (\beta B' P A + N) ``` The sequence $\{d_t\}$ from {eq}`lq_dd` is replaced by the constant value @@ -904,7 +899,7 @@ The sequence $\{d_t\}$ from {eq}`lq_dd` is replaced by the constant value :label: lq_dd_ih d -:= \mathop{\mathrm{trace}}(C^\top P C) \frac{\beta}{1 - \beta} +:= \mathop{\mathrm{trace}}(C' P C) \frac{\beta}{1 - \beta} ``` The state evolves according to the time-homogeneous process $x_{t+1} = (A - BF) x_t + C w_{t+1}$. @@ -977,36 +972,44 @@ Thus, for the dynamics we set :label: lq_lowmc3 x_t := -\begin{bmatrix} +\left( +\begin{array}{c} a_t \\ 1 \\ t \\ t^2 -\end{bmatrix}, +\end{array} +\right), \quad A := -\begin{bmatrix} +\left( +\begin{array}{cccc} 1 + r & -\bar c & m_1 & m_2 \\ 0 & 1 & 0 & 0 \\ 0 & 1 & 1 & 0 \\ 0 & 1 & 2 & 1 -\end{bmatrix}, +\end{array} +\right), \quad B := -\begin{bmatrix} +\left( +\begin{array}{c} -1 \\ 0 \\ 0 \\ 0 -\end{bmatrix}, +\end{array} +\right), \quad C := -\begin{bmatrix} +\left( +\begin{array}{c} \sigma \\ 0 \\ 0 \\ 0 -\end{bmatrix} +\end{array} +\right) ``` If you expand the expression $x_{t+1} = A x_t + B u_t + C w_{t+1}$ using @@ -1021,20 +1024,24 @@ To implement preference specification {eq}`lq_pip` we take Q := 1, \quad R := -\begin{bmatrix} +\left( +\begin{array}{cccc} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 -\end{bmatrix} +\end{array} +\right) \quad \text{and} \quad R_f := -\begin{bmatrix} +\left( +\begin{array}{cccc} q & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 -\end{bmatrix} +\end{array} +\right) ``` The next figure shows a simulation of consumption and assets computed using @@ -1203,7 +1210,7 @@ To produce these figures we converted the monopolist problem into an LQ problem. The key to this conversion is to choose the right state --- which can be a bit of an art. -Here we take $x_t = (\bar q_t \;\, q_t \;\, 1)^\top$, while the control is chosen as $u_t = q_{t+1} - q_t$. +Here we take $x_t = (\bar q_t \;\, q_t \;\, 1)'$, while the control is chosen as $u_t = q_{t+1} - q_t$. We also manipulated the profit function slightly. @@ -1256,7 +1263,6 @@ The parameters are $r = 0.05, \beta = 1 / (1 + r), \bar c = 1.5, \mu = 2, \sigm ```{solution-start} lqc_ex1 :class: dropdown -:label: lqc_ex1_solution ``` Here’s one solution. @@ -1386,7 +1392,6 @@ together the simulations from these two separate models. ```{solution-start} lqc_ex2 :class: dropdown -:label: lqc_ex2_solution ``` This is a permanent income / life-cycle model with polynomial growth in @@ -1510,13 +1515,12 @@ For parameters, use $a_0 = 5, a_1 = 0.5, \sigma = 0.15, \rho = 0.9, ```{solution-start} lqc_ex3 :class: dropdown -:label: lqc_ex3_solution ``` The first task is to find the matrices $A, B, C, Q, R$ that define the LQ problem. -Recall that $x_t = (\bar q_t \;\, q_t \;\, 1)^\top$, while +Recall that $x_t = (\bar q_t \;\, q_t \;\, 1)'$, while $u_t = q_{t+1} - q_t$. Letting $m_0 := (a_0 - c) / 2a_1$ and $m_1 := 1 / 2 a_1$, we