From 15b8ebb8612d94a84b87f1f5549b6ec6aee27014 Mon Sep 17 00:00:00 2001 From: thomassargent30 Date: Fri, 31 Jul 2026 16:52:30 -0600 Subject: [PATCH 1/3] Tom's July 31 edits of some old and new lectures --- lectures/_static/quant-econ.bib | 225 +++ lectures/_toc.yml | 5 +- lectures/ar1_turningpts.md | 2 +- lectures/kalman_filter_var.md | 434 +---- lectures/lq_bewley_complete_markets.md | 6 +- lectures/lq_permanent_income.md | 5 +- lectures/lq_robust_bewley.md | 622 +++++++ lectures/lq_robust_smoothing.md | 1272 +++++++------- lectures/sargent_surico.md | 2096 ++++++++++++++++++++++++ lectures/var_subsets.md | 1234 ++++++++++++++ 10 files changed, 4938 insertions(+), 963 deletions(-) create mode 100644 lectures/lq_robust_bewley.md create mode 100644 lectures/sargent_surico.md create mode 100644 lectures/var_subsets.md diff --git a/lectures/_static/quant-econ.bib b/lectures/_static/quant-econ.bib index 2b65618fb..9594beaee 100644 --- a/lectures/_static/quant-econ.bib +++ b/lectures/_static/quant-econ.bib @@ -5000,3 +5000,228 @@ @incollection{Axelrod1987 pages = {32--41}, year = {1987} } + +@article{SargentSurico2011, + author = {Sargent, Thomas J. and Surico, Paolo}, + title = {Two Illustrations of the Quantity Theory of Money: Breakdowns and Revivals}, + journal = {American Economic Review}, + volume = {101}, + number = {1}, + pages = {109--128}, + year = {2011}, + doi = {10.1257/aer.101.1.109} +} + +@article{Lucas1980, + author = {Lucas, Robert E., Jr.}, + title = {Two Illustrations of the Quantity Theory of Money}, + journal = {American Economic Review}, + volume = {70}, + number = {5}, + pages = {1005--1014}, + year = {1980} +} + +@article{Whiteman1984, + author = {Whiteman, Charles H.}, + title = {Lucas on the Quantity Theory: Hypothesis Testing without Theory}, + journal = {American Economic Review}, + volume = {74}, + number = {4}, + pages = {742--749}, + year = {1984} +} + +@article{Ireland2004, + author = {Ireland, Peter N.}, + title = {Technology Shocks in the New Keynesian Model}, + journal = {Review of Economics and Statistics}, + volume = {86}, + number = {4}, + pages = {923--936}, + year = {2004} +} + +@article{Ireland2003, + author = {Ireland, Peter N.}, + title = {Endogenous Money or Sticky Prices?}, + journal = {Journal of Monetary Economics}, + volume = {50}, + number = {8}, + pages = {1623--1648}, + year = {2003} +} + +@article{Sims2002gensys, + author = {Sims, Christopher A.}, + title = {Solving Linear Rational Expectations Models}, + journal = {Computational Economics}, + volume = {20}, + number = {1--2}, + pages = {1--20}, + year = {2002} +} + +@article{AnSchorfheide2007, + author = {An, Sungbae and Schorfheide, Frank}, + title = {Bayesian Analysis of {DSGE} Models}, + journal = {Econometric Reviews}, + volume = {26}, + number = {2--4}, + pages = {113--172}, + year = {2007} +} + +@article{LubikSchorfheide2004, + author = {Lubik, Thomas A. and Schorfheide, Frank}, + title = {Testing for Indeterminacy: An Application to {U.S.} Monetary Policy}, + journal = {American Economic Review}, + volume = {94}, + number = {1}, + pages = {190--217}, + year = {2004} +} + +@incollection{McCallumNelson1999, + author = {McCallum, Bennett T. and Nelson, Edward}, + title = {Performance of Operational Policy Rules in an Estimated Semiclassical Structural Model}, + editor = {Taylor, John B.}, + booktitle = {Monetary Policy Rules}, + publisher = {University of Chicago Press}, + address = {Chicago}, + pages = {15--45}, + year = {1999} +} + +@article{Rotemberg1982, + author = {Rotemberg, Julio J.}, + title = {Sticky Prices in the United States}, + journal = {Journal of Political Economy}, + volume = {90}, + number = {6}, + pages = {1187--1211}, + year = {1982} +} + +@incollection{BalkeGordon1986, + author = {Balke, Nathan S. and Gordon, Robert J.}, + title = {Appendix B: Historical Data}, + editor = {Gordon, Robert J.}, + booktitle = {The American Business Cycle: Continuity and Change}, + publisher = {University of Chicago Press}, + address = {Chicago}, + pages = {781--850}, + year = {1986} +} + +@book{FriedmanSchwartz1963, + author = {Friedman, Milton and Schwartz, Anna J.}, + title = {A Monetary History of the United States, 1867--1960}, + publisher = {Princeton University Press}, + address = {Princeton, NJ}, + year = {1963} +} + +@article{SmetsWouters2007, + author = {Smets, Frank and Wouters, Rafael}, + title = {Shocks and Frictions in {US} Business Cycles: A {Bayesian} {DSGE} Approach}, + journal = {American Economic Review}, + volume = {97}, + number = {3}, + pages = {586--606}, + year = {2007} +} + +@article{Sargent1971, + author = {Sargent, Thomas J.}, + title = {A Note on the `Accelerationist' Controversy}, + journal = {Journal of Money, Credit and Banking}, + volume = {3}, + number = {3}, + pages = {721--725}, + year = {1971} +} + +@article{Lucas1975, + author = {Lucas, Robert E., Jr.}, + title = {An Equilibrium Model of the Business Cycle}, + journal = {Journal of Political Economy}, + volume = {83}, + number = {6}, + pages = {1113--1144}, + year = {1975} +} + +@article{BoschenOtrok1994, + author = {Boschen, John F. and Otrok, Christopher M.}, + title = {Long-Run Neutrality and Superneutrality in an {ARIMA} Framework: Comment}, + journal = {American Economic Review}, + volume = {84}, + number = {5}, + pages = {1470--1473}, + year = {1994} +} + +@techreport{SargentSurico2008, + author = {Sargent, Thomas J. and Surico, Paolo}, + title = {Monetary Policies and Low-Frequency Manifestations of the Quantity Theory}, + institution = {Bank of England External Monetary Policy Committee Unit}, + type = {Discussion Paper}, + number = {26}, + year = {2008} +} + +@article{DuaneEtAl1987, + author = {Duane, Simon and Kennedy, A. D. and Pendleton, Brian J. and Roweth, Duncan}, + title = {Hybrid {Monte} {Carlo}}, + journal = {Physics Letters B}, + volume = {195}, + number = {2}, + pages = {216--222}, + year = {1987}, + doi = {10.1016/0370-2693(87)91197-X} +} + +@incollection{Neal2011, + author = {Neal, Radford M.}, + title = {{MCMC} Using {Hamiltonian} Dynamics}, + editor = {Brooks, Steve and Gelman, Andrew and Jones, Galin L. and Meng, Xiao-Li}, + booktitle = {Handbook of {Markov} Chain {Monte} {Carlo}}, + publisher = {Chapman and Hall/CRC}, + pages = {113--162}, + year = {2011} +} + +@article{HoffmanGelman2014, + author = {Hoffman, Matthew D. and Gelman, Andrew}, + title = {The No-U-Turn Sampler: Adaptively Setting Path Lengths in {Hamiltonian} {Monte} {Carlo}}, + journal = {Journal of Machine Learning Research}, + volume = {15}, + number = {47}, + pages = {1593--1623}, + year = {2014} +} + +@article{Betancourt2017, + author = {Betancourt, Michael}, + title = {A Conceptual Introduction to {Hamiltonian} {Monte} {Carlo}}, + journal = {arXiv preprint arXiv:1701.02434}, + year = {2017} +} + +@article{Klein2000, + author = {Klein, Paul}, + title = {Using the Generalized {Schur} Form to Solve a Multivariate Linear Rational Expectations Model}, + journal = {Journal of Economic Dynamics and Control}, + volume = {24}, + number = {10}, + pages = {1405--1423}, + year = {2000} +} + +@article{PhanEtAl2019, + author = {Phan, Du and Pradhan, Neeraj and Jankowiak, Martin}, + title = {Composable Effects for Flexible and Accelerated Probabilistic Programming in {NumPyro}}, + journal = {arXiv preprint arXiv:1912.11554}, + year = {2019} +} diff --git a/lectures/_toc.yml b/lectures/_toc.yml index 91397a3c1..dcd0a4545 100644 --- a/lectures/_toc.yml +++ b/lectures/_toc.yml @@ -65,7 +65,8 @@ parts: - file: wealth_dynamics - file: kalman - file: kalman_2 - - file: kalman_filter_var + - file: kalman_filter_var + - file: var_subsets - file: organization_capital - file: measurement_models - caption: Search @@ -116,6 +117,7 @@ parts: - file: lq_permanent_income - file: lq_bewley_complete_markets - file: lq_robust_smoothing + - file: lq_robust_bewley - file: lq_inventories - caption: Bounded Rationality in Macroeconomics numbered: true @@ -185,6 +187,7 @@ parts: - file: mle - file: unemployment_linear - file: unemployment_shocks + - file: sargent_surico - caption: Auctions numbered: true chapters: diff --git a/lectures/ar1_turningpts.md b/lectures/ar1_turningpts.md index 1191d15e0..b21e015d4 100644 --- a/lectures/ar1_turningpts.md +++ b/lectures/ar1_turningpts.md @@ -505,7 +505,7 @@ def draw_from_posterior(data, size=10000, dis_plot=True, key=key): # Plot posterior distributions and trace plots if dis_plot: plot_data = az.from_numpyro(posterior=mcmc) - az.plot_trace_dist(plot_data, var_names=['ρ', 'σ']) + az.plot_trace(plot_data, var_names=['ρ', 'σ']) return post_sample diff --git a/lectures/kalman_filter_var.md b/lectures/kalman_filter_var.md index b0e0c79ef..52a61a773 100644 --- a/lectures/kalman_filter_var.md +++ b/lectures/kalman_filter_var.md @@ -4,13 +4,22 @@ jupytext: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.11.1 + jupytext_version: 1.16.7 kernelspec: - display_name: Python 3 + display_name: Python 3 (ipykernel) language: python name: python3 --- +(kalman_filter_var)= +```{raw} jupyter +
+ + QuantEcon + +
+``` + # The Kalman Filter and Vector Autoregressions ```{index} single: Kalman Filter @@ -19,6 +28,10 @@ kernelspec: ```{index} single: Vector Autoregression; and Kalman filter ``` +```{contents} Contents +:depth: 2 +``` + In addition to what's in Anaconda, this lecture will need the following libraries: ```{code-cell} ipython3 @@ -54,6 +67,10 @@ The lecture covers: - why the Kalman filter is an essential tool for *interpreting VARs* estimated from economic data +A sequel, {doc}`var_subsets`, applies this machinery to a question that arises +constantly in practice: what happens to a VAR when an econometrician observes only +some of the variables that appear in it. + ## The state space system The Kalman filter applies to the **state space system** for $t \geq 0$: @@ -73,7 +90,7 @@ where $0$ and identity covariance matrix - $v_t$ is an IID sequence of normal random variables with mean zero and covariance matrix $R$ -- $w_{t+1}$ and $v_s$ are orthogonal for all $t+1$ and $s \geq 0$ +- $w_t$ and $v_s$ are orthogonal at all pairs of dates The coefficient matrices have the following dimensions: $A$ is $n \times n$, $C$ is $n \times p$, $G$ is $m \times n$, and $R$ is $m \times m$. @@ -107,7 +124,7 @@ $y_t$ given history $y^{t-1}$. The Kalman filter attains this by constructing recursive formulas for $\hat{x}_t$ and $\Sigma_t$ such that the distribution of $y_t$ conditional on -$y^{t-1}$ generalises {eq}`eq:kalf4` to +$y^{t-1}$ generalizes {eq}`eq:kalf4` to $$ y_t \sim N(G \hat{x}_t,\; G \Sigma_t G^\top + R) @@ -116,7 +133,7 @@ $$ (eq:kalf400) for $t \geq 1$, where the distribution of $x_t$ conditional on $y^{t-1}$ is $N(\hat{x}_t, \Sigma_t)$. -The objects $\hat{x}_t$ and $\Sigma_t$ characterise the **population regression** +The objects $\hat{x}_t$ and $\Sigma_t$ characterize the **population regression** $$ \hat{x}_t = \mathbb{E}[x_t \mid y_{t-1}, \ldots, y_0] @@ -231,7 +248,7 @@ $$ (eq:kalf1000) System {eq}`eq:kalf1000` maps a mean-covariance pair $(\hat{x}_0, \Sigma_0)$ into a new pair $(\hat{x}_1, \Sigma_1)$, with auxiliary outputs $(a_0, K_0)$. -Recognising that "we are in the same situation at the start of period 1 as at +Recognizing that "we are in the same situation at the start of period 1 as at the start of period 0" activates a recursion, the **Kalman filter**. ### The Kalman filter recursions @@ -297,8 +314,8 @@ Thus $\{a_t\}$ is a white-noise process of innovations to $\{y_t\}$. Sometimes {eq}`eq:kalf10` is called a **whitening filter**: it takes the signal process $\{y_t\}$ as input and produces the white-noise innovation process $\{a_t\}$ as output. -With $H(a^t)$ defined analogously, the linear space $H(a^t)$ is an orthogonal -basis for the linear space $H(y^t)$. +With $H(a^t)$ defined analogously, $H(a^t) = H(y^t)$, and $[a_t, \ldots, a_0]$ +is an orthogonal basis for that common space. Rather than computing $\mathbb{E}[x_t \mid y_{t-1}, \ldots, y_0]$ via one large regression, the Kalman filter performs a sequence of small regressions on successive orthogonal @@ -342,7 +359,7 @@ For $t \geq 1$, $\mathbb{E}[y_t \mid y^{t-1}] = G\hat{x}_t$ and the conditional distribution of $y_t$ given $y^{t-1}$ is $N(G\hat{x}_t, \Omega_t)$. The objects $(G\hat{x}_t, \Omega_t)$ emerging from the Kalman filter recursions -therefore completely characterise this conditional distribution. +therefore completely characterize this conditional distribution. ### The likelihood function @@ -405,7 +422,7 @@ where the conditioning extends over the **semi-infinite** past $s \leq t-1$. ### A time-invariant VAR -If the fixed point $\Sigma$ exists and we initialise the filter at $\Sigma_0 = \Sigma$, +If the fixed point $\Sigma$ exists and we initialize the filter at $\Sigma_0 = \Sigma$, the innovations representation {eq}`eq:innovrep` becomes time-invariant: $$ @@ -443,7 +460,7 @@ $$ (eq:varorth) The orthogonality conditions {eq}`eq:varorth` identify {eq}`eq:var1` as a vector autoregression. -Defining the lag operator $L$ by $L x_{t+1} \equiv x_t$, the +Letting $L$ denote the lag operator, so that $L x_t = x_{t-1}$, the **moving average representation** deduced from {eq}`eq:innovti` is $$ @@ -647,24 +664,37 @@ We first simulate a sample path of the true hidden state and noisy observations. T = 200 x_path, y_path = lss.simulate(ts_length=T, random_state=42) -# Shapes: x_path is (n, T+1), y_path is (m, T) -x_true = x_path[0, :T] +# Shapes: x_path is (n, T), y_path is (m, T) +x_true = x_path[0, :] y_obs = y_path[0, :] ``` -We then run the Kalman filter manually, step by step, to collect the filtered -estimates. +We then run the Kalman filter manually, step by step. + +The `Kalman.update` method performs a *complete* cycle, moving the prior to the +filtering distribution and then on to next period's prior. + +So $(\hat{x}_t, \Sigma_t)$ as defined in {eq}`eq:kalf10` are the values held by +the object *before* `update` is called, and we record them there. ```{code-cell} ipython3 x_hats = np.zeros(T) Sigmas = np.zeros(T) +innovations = np.zeros(T) for t in range(T): + x_hats[t] = kf.x_hat.item() # x_hat_t = E[x_t | y^{t-1}] + Sigmas[t] = kf.Sigma.item() # Sigma_t + innovations[t] = y_obs[t] - (G @ kf.x_hat).item() kf.update(y_obs[t:t+1]) # one full filter cycle - x_hats[t] = kf.x_hat.item() - Sigmas[t] = kf.Sigma.item() ``` +Getting this ordering right matters. + +Recording `kf.x_hat` *after* the call would store the one-step-ahead forecast +$\hat{x}_{t+1}$, and differencing it from $y_t$ would produce a series that is +not the innovation and does not have variance $G\Sigma G^\top + R$. + ```{code-cell} ipython3 --- mystnb: @@ -691,7 +721,7 @@ axes[1].axhline(kf.Sigma_infinity[0, 0], ls='--', color='k', axes[1].set_title('conditional variance') axes[1].legend(fontsize=9) -axes[2].plot(t_range, y_obs - x_hats, color='C2', lw=2, alpha=0.7, +axes[2].plot(t_range, innovations, color='C2', lw=2, alpha=0.7, label=r'innovation $a_t = y_t - G\hat{x}_t$') axes[2].set_title('innovation') axes[2].set_xlabel('time $t$') @@ -709,6 +739,19 @@ onto its steady-state value. The innovation series fluctuates around zero, as a one-step forecast error should. +Its sample standard deviation should be close to $\sqrt{G\Sigma_\infty G^\top + R}$. + +```{code-cell} ipython3 +print(f"sample sd of innovations = {innovations.std():.4f}") +print(f"steady-state sd = " + f"{np.sqrt(kf.Sigma_infinity[0, 0] + R[0, 0]):.4f}") +print(f"first-order autocorrelation = " + f"{np.corrcoef(innovations[1:], innovations[:-1])[0, 1]:.4f}") +``` + +The autocorrelation is near zero, confirming that the filter has whitened the +observed series. + ### Convergence of the Riccati equation The `Kalman` class computes the steady-state covariance $\Sigma_\infty$ by @@ -797,342 +840,23 @@ ll = log_likelihood(A, C, G, R, print(f"Log-likelihood of sample: {ll:.4f}") ``` -## An example - -We now work through a structured example that shows how a bivariate VAR(2) fits naturally into the state space framework and how the Kalman filter delivers a -Wold (innovations) representation. - -### A linear state-space system and its filter - -The state and observation equations are - -$$ -x_{t+1} = A x_t + C w_{t+1} -$$ (eq:ex_state) - -$$ -y_t = G x_t + v_t -$$ (eq:ex_obs) - -with initial condition and shock distributions - -$$ -x_0 \sim N(\hat{x}_0, \Sigma_0), \quad -w_{t+1} \sim N(0, I), \quad -v_t \sim N(0, R). -$$ - -The steady-state error covariance matrix $\Sigma$ satisfies the Riccati equation - -$$ -\Sigma = A \Sigma A^\top + CC^\top - - A \Sigma G^\top \bigl(G \Sigma G^\top + R\bigr)^{-1} G \Sigma A^\top -$$ (eq:ex_riccati) - -and the associated steady-state Kalman gain is +## Where this leads -$$ -K = A \Sigma G^\top \bigl(G \Sigma G^\top + R\bigr)^{-1} -$$ (eq:ex_gain) - -Starting from an initial estimate $\hat{x}_0$, the Kalman filter updates the state estimate via - -$$ -\hat{x}_{t+1} = A \hat{x}_t + K a_t -$$ (eq:ex_kf_update) - -where the innovation is - -$$ -a_t = y_t - G \hat{x}_t -$$ (eq:ex_innovation) - -Substituting {eq}`eq:ex_innovation` into {eq}`eq:ex_kf_update` and expanding: +The Kalman filter maps a state space system into the VAR that a population +regression of $y_t$ on its own past would recover. -$$ -\hat{x}_{t+1} = A \hat{x}_t + K(y_t - G\hat{x}_t) - = (A - KG)\hat{x}_t + K y_t - = (A - KG)\hat{x}_t + K G x_t + K v_t -$$ (eq:ex_kf_expanded) - -### Impulse responses of $y_t$ to the innovations $a_t$ - -It is useful to compute the **ordinary impulse response functions** of the -observable vector $y_t$ to its own innovations $a_t$, the moving-average (Wold) -representation that is the mirror image of the VAR {eq}`eq:var1`. - -From the time-invariant innovations representation {eq}`eq:innovti` - -$$ -\hat{x}_{t+1} = A\hat{x}_t + K a_t, \qquad y_t = G\hat{x}_t + a_t, -$$ +That map is most interesting when the state $x_t$ contains variables that the +econometrician does *not* see. -the moving-average representation {eq}`eq:sf_wold` is - -$$ -y_t = \bigl[I + G(I - AL)^{-1} K L\bigr]\, a_t - = a_t + \sum_{h=1}^{\infty} G A^{h-1} K\, a_{t-h}. -$$ - -Hence the impulse response of $y_t$ to a unit innovation $a_t$ is - -$$ -\Psi_0 = I, \qquad \Psi_h = G A^{h-1} K \quad (h \ge 1). -$$ (eq:ex_y_to_a) - -These coefficients decay at the rate governed by the eigenvalues of $A$. - -We can read the coefficients {eq}`eq:ex_y_to_a` directly off a `quantecon` -`LinearStateSpace` object. - -We build a state-space system whose state is the -filtered estimate $\hat{x}_t$, whose single "shock" is the innovation $a_t$ -loaded through $C = K$, and whose observation matrix is $G$. - -The -`impulse_response` method of that object returns the sequence $G A^{j} K$ for -$j = 0, 1, 2, \ldots$, which are exactly the $\Psi_h$ for $h \ge 1$; we prepend -$\Psi_0 = I$ to capture the contemporaneous feed-through $y_t = G\hat{x}_t + a_t$. - -The array returned below has entry `[h, i, j]` equal to the response of -observable `i` at horizon `h` to innovation component `j`. - -```{code-cell} ipython3 -def y_to_a_irf(A, K, G, T=40): - """ - Return Wold IRFs of y_t to its own innovations a_t. - """ - n, m = A.shape[0], G.shape[0] - lss = qe.LinearStateSpace(A, K, G, np.zeros((m, m)), mu_0=np.zeros(n)) - _, ycoef = lss.impulse_response(j=T - 2) # [GK, GAK, GA^2K, ...] - Psi = np.empty((T, m, m)) - Psi[0] = np.eye(m) # contemporaneous response - for h in range(1, T): - Psi[h] = ycoef[h - 1] - return Psi -``` +The sequel {doc}`var_subsets` takes the leading case of this: $Y_t$ follows a +finite-order VAR, and the econometrician observes only a subvector +$y_t = S_y Y_t$. -### Bivariate VAR(2) in state-space form +There we build general code that, given the VAR for $Y_t$ and the selector +matrix $S_y$, returns the VAR and moving average representations for $y_t$ and +expresses the innovations in the small system as a distributed lag of the +innovations in the large one. -Consider two observable series $r_t$ and $z_t$. - -Stack them into the state vector $x_t = (r_t,\; r_{t-1},\; z_t,\; z_{t-1})^\top$. - -We posit the VAR(2) state-transition equation: - -$$ -\begin{pmatrix} r_{t+1} \\ r_t \\ z_{t+1} \\ z_t \end{pmatrix} -= -\begin{pmatrix} - d_1 & d_2 & d_3 & d_4 \\ - 1 & 0 & 0 & 0 \\ - \delta_1 & \delta_2 & \delta_3 & \delta_4 \\ - 0 & 0 & 1 & 0 -\end{pmatrix} -\begin{pmatrix} r_t \\ r_{t-1} \\ z_t \\ z_{t-1} \end{pmatrix} -+ -\begin{pmatrix} - c_{11} & c_{12} \\ - 0 & 0 \\ - c_{21} & c_{22} \\ - 0 & 0 -\end{pmatrix} -\begin{pmatrix} w_{1,t+1} \\ w_{2,t+1} \end{pmatrix} -$$ (eq:ex_var2_state) - -We consider two possible observation equations. - -The first is a bivariate observation of $r_t$ and $z_t$: - -$$ -\begin{pmatrix} r_t \\ z_t \end{pmatrix} -= -\begin{pmatrix} - 1 & 0 & 0 & 0 \\ - 0 & 0 & 1 & 0 -\end{pmatrix} -\begin{pmatrix} r_t \\ r_{t-1} \\ z_t \\ z_{t-1} \end{pmatrix} -+ -\begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} -\begin{pmatrix} v_{1t} \\ v_{2t} \end{pmatrix} -$$ (eq:ex_var2_obs) - -The second is a univariate observation of $r_t$: - -$$ -y_t = \begin{pmatrix} 1 & 0 & 0 & 0 \end{pmatrix} -\begin{pmatrix} r_t \\ r_{t-1} \\ z_t \\ z_{t-1} \end{pmatrix} -+ v_{1t} -$$ (eq:ex_scalar_obs) - -We now compare the Wold impulse responses generated by these two observation -systems. - -System 1 observes both $r_t$ and $z_t$, so its innovation $a_t$ is -$2 \times 1$. - -System 2 observes only $r_t$, so its innovation $u_t$ is scalar. - -The transition matrices are the same in the two systems, but the observation -matrix changes. - -Consequently, the steady-state Kalman gain changes too, and so do the Wold -responses of the observables to their own innovations. - -### Numerical example: impulse responses to innovations - -The parameter values are: - -$$ -\begin{aligned} -d_1 &= 0.80,\quad d_2 = 0.05,\quad d_3 = 0.75,\quad d_4 = -0.72 \\ -\delta_1 &= 0.00,\quad \delta_2 = 0.00,\quad \delta_3 = 0.75,\quad \delta_4 = 0.20 \\ -c_{11} &= 1.0,\quad c_{12} = 0.0,\quad c_{21} = 0.0,\quad c_{22} = 1.0 \\ -R &= 0.0001 \times I_2 \quad \text{(bivariate case)}, \qquad -R = 0.0001 \quad \text{(univariate case)}. -\end{aligned} -$$ - -These give the $4 \times 4$ transition matrix and $4 \times 2$ shock-loading matrix - -$$ -A = \begin{pmatrix} -0.80 & 0.05 & 0.75 & -0.72 \\ -1 & 0 & 0 & 0 \\ -0 & 0 & 0.75 & 0.20 \\ -0 & 0 & 1 & 0 -\end{pmatrix}, \qquad -C = \begin{pmatrix} -1 & 0 \\ -0 & 0 \\ -0 & 1 \\ -0 & 0 -\end{pmatrix}. -$$ - -**System 1** uses the bivariate observation equation {eq}`eq:ex_var2_obs`, so -$G$ selects $(r_t, z_t)^\top$ from the state and the innovation $a_t$ is $2 \times 1$. - -**System 2** uses the univariate observation equation {eq}`eq:ex_scalar_obs`, so -the row vector in that equation selects only $r_t$ and the innovation $u_t$ is scalar. - -```{code-cell} ipython3 -# Parameters -d1, d2, d3, d4 = 0.80, 0.05, 0.75, -.72 -δ1, δ2, δ3, δ4 = 0.00, 0.00, 0.75, 0.20 -c11, c12, c21, c22 = 1.0, 0.0, 0.0, 1.0 -σ_v = 0.01 # sqrt(0.0001) - -# Shared matrices -A_var = np.array([[d1, d2, d3, d4 ], - [1.0, 0.0, 0.0, 0.0 ], - [δ1, δ2, δ3, δ4], - [0.0, 0.0, 1.0, 0.0 ]]) - -C_var = np.array([[c11, c12], - [0.0, 0.0], - [c21, c22], - [0.0, 0.0]]) - -# System 1: bivariate observation -G_biv = np.array([[1.0, 0.0, 0.0, 0.0], - [0.0, 0.0, 1.0, 0.0]]) -H_biv = σ_v * np.eye(2) # H @ H.T = 0.0001 * I_2 - -lss_biv = qe.LinearStateSpace(A_var, C_var, G_biv, H_biv, - mu_0=np.zeros(4), Sigma_0=np.eye(4)) -kf_biv = qe.Kalman(lss_biv) -_, K_biv = kf_biv.stationary_values() - -print("System 1 - steady-state Kalman gain K (4x2):") -print(np.round(K_biv, 5)) - -# System 2: univariate observation -G_uni = np.array([[1.0, 0.0, 0.0, 0.0]]) -H_uni = np.array([[σ_v]]) # H @ H.T = 0.0001 - -lss_uni = qe.LinearStateSpace(A_var, C_var, G_uni, H_uni, - mu_0=np.zeros(4), Sigma_0=np.eye(4)) -kf_uni = qe.Kalman(lss_uni) -_, K_uni = kf_uni.stationary_values() - -print("\nSystem 2 - steady-state Kalman gain K (4x1):") -print(np.round(K_uni, 5)) -``` - -We now apply the helper `y_to_a_irf` defined above to compute the ordinary -impulse responses {eq}`eq:ex_y_to_a` of the observable $y_t$ to its own -innovations $a_t$, for both System 1 (bivariate, so $a_t$ is $2 \times 1$) and -System 2 (univariate, so $u_t$ is scalar). - -```{code-cell} ipython3 ---- -mystnb: - figure: - caption: System 1 responses to own innovations - name: fig-kfvar-sys1-ya ---- -T_irf = 40 -horizons = np.arange(T_irf) - -Psi_biv = y_to_a_irf(A_var, K_biv, G_biv, T_irf) # System 1: (T, 2, 2) -Psi_uni = y_to_a_irf(A_var, K_uni, G_uni, T_irf) # System 2: (T, 1, 1) - -obs_labels = [r'$r_t$', r'$z_t$'] -innov_labels = [r'$a_{1,t}$', r'$a_{2,t}$'] - -# System 1 responses -fig, axes = plt.subplots(2, 2, figsize=(10, 6), sharex=True) -for i, obs in enumerate(obs_labels): - for j, inn in enumerate(innov_labels): - ax = axes[i, j] - ax.plot(horizons, Psi_biv[:, i, j], lw=2) - ax.axhline(0, color='k', lw=0.6, ls='--') - ax.set_title(fr'{obs} to {inn}', fontsize=9) - if i == 1: - ax.set_xlabel('horizon $h$') - if j == 0: - ax.set_ylabel('response') -fig.tight_layout() -plt.show() -``` - -Own innovations move their own observables one-for-one on impact and then fade. - -The diagonal panels start at one and the off-diagonal panels start at zero -because $\Psi_0 = I$. - -The cross response of $r_t$ to $a_{2,t}$ is sizeable at short horizons, while -the response of $z_t$ to $a_{1,t}$ is tiny on the displayed scale. - -```{code-cell} ipython3 ---- -mystnb: - figure: - caption: System 2 response to own innovation - name: fig-kfvar-sys2-ya ---- -fig, ax = plt.subplots() -ax.plot(horizons, Psi_uni[:, 0, 0], lw=2) -ax.axhline(0, color='k', lw=0.6, ls='--') -ax.set_xlabel('horizon $h$') -ax.set_ylabel('response') -fig.tight_layout() -plt.show() -``` - -With only $r_t$ observed, the single innovation moves $r_t$ one-for-one on -impact and then decays monotonically toward zero. - -For $h \ge 1$ the responses propagate through the state matrix $A$ and decay -geometrically, tracing out the Wold moving-average representation of the -bivariate (System 1) and univariate (System 2) processes. - -These are forecast-error responses from the Wold representation, not structural -shock responses. - -With the example in place, we can now pull together the main lessons of the -lecture. ## Summary @@ -1149,9 +873,11 @@ construction. Solving the innovations representation backward gives an infinite-order VAR, while solving it forward gives the Wold moving-average representation. -The numerical examples show that changing the observed variables changes the -Kalman gain and therefore changes the Wold responses, even when the underlying -state dynamics are the same. +Which variables an econometrician observes determines the Kalman gain and +therefore the VAR and Wold representations, even when the underlying state +dynamics are held fixed. + +{doc}`var_subsets` pursues that observation systematically. ## Exercises @@ -1278,12 +1004,12 @@ x2_path, y2_path = lss2.simulate(ts_length=T2, random_state=0) x_hats2 = np.zeros((T2, 2)) for t in range(T2): + x_hats2[t] = kf2.x_hat.ravel() # record x_hat_t before updating kf2.update(y2_path[:, t]) - x_hats2[t] = kf2.x_hat.ravel() fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True) for i, ax in enumerate(axes): - ax.plot(x2_path[i, :T2], lw=2, label=f'true $x_{{{i+1},t}}$') + ax.plot(x2_path[i, :], lw=2, label=f'true $x_{{{i+1},t}}$') ax.plot(x_hats2[:, i], lw=2, ls='--', label=rf'$\hat{{x}}_{{{i+1},t}}$') ax.set_title(f'component {i+1}') ax.legend(fontsize=9) diff --git a/lectures/lq_bewley_complete_markets.md b/lectures/lq_bewley_complete_markets.md index 70352ddce..ba867623f 100644 --- a/lectures/lq_bewley_complete_markets.md +++ b/lectures/lq_bewley_complete_markets.md @@ -39,7 +39,7 @@ kernelspec: This lecture studies how the cross-section distribution of consumption evolves when many consumers each solve the LQ permanent income problem. -It is the second of three lectures on the LQ permanent income model and builds directly on {doc}`lq_permanent_income`. +It is the second of four lectures on the LQ permanent income model and builds directly on {doc}`lq_permanent_income`. We first show that the unit root in individual consumption causes the cross-section variance of consumption to grow linearly with time. @@ -49,6 +49,8 @@ Finally, we replace the single risk-free bond with a complete set of Arrow secur The third lecture, {doc}`lq_robust_smoothing`, relaxes the assumption that the consumer fully trusts his income model. +The fourth, {doc}`lq_robust_bewley`, returns to the Bewley economy of this lecture and populates it with consumers who differ in how much they trust that model. + Let's begin with some imports. ```{code-cell} ipython3 @@ -439,6 +441,8 @@ So far the consumer fully trusts his stochastic income model. In {doc}`lq_robust_smoothing` we relax that assumption and let the consumer seek decision rules that are robust to plausible misspecifications. +{doc}`lq_robust_bewley` then returns to the Bewley economy studied here and shows that a cross-section of consumers with different concerns about misspecification can reproduce exactly the equilibrium of this lecture. + The optimal robust rule takes the same form as the rule above, but under a distorted model of the income process that looks more persistent than the approximating one. ## Exercises diff --git a/lectures/lq_permanent_income.md b/lectures/lq_permanent_income.md index 4afbb9f28..5e1d139fa 100644 --- a/lectures/lq_permanent_income.md +++ b/lectures/lq_permanent_income.md @@ -48,12 +48,13 @@ The model is useful for studying We derive the consumer's optimal consumption function, present two state-space representations of the optimal decision rule, and illustrate them with two classic examples. -This is the first of three lectures on the LQ permanent income model. +This is the first of four lectures on the LQ permanent income model. The two sequels build directly on the tools developed here. - {doc}`lq_bewley_complete_markets` studies the cross-section behavior of consumption and embeds the single consumer in closed economies with incomplete and complete markets. - {doc}`lq_robust_smoothing` studies a consumer who distrusts his model of income and engages in precautionary savings. +- {doc}`lq_robust_bewley` combines the two, building a Bewley economy whose consumers differ in how much they distrust their income model. Let's begin with some imports. @@ -510,6 +511,8 @@ The two representations and examples developed here are the foundation for the t {doc}`lq_robust_smoothing` studies a consumer who distrusts the endowment process {eq}`eq:sprob15` and engages in precautionary savings. +{doc}`lq_robust_bewley` puts such consumers into a Bewley economy and shows that their concerns about misspecification leave no trace in quantity data. + ## Exercises ```{exercise-start} diff --git a/lectures/lq_robust_bewley.md b/lectures/lq_robust_bewley.md new file mode 100644 index 000000000..212ea86e3 --- /dev/null +++ b/lectures/lq_robust_bewley.md @@ -0,0 +1,622 @@ +--- +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.11.1 +kernelspec: + display_name: Python 3 + language: python + name: python3 +--- + +(lq_robust_bewley)= +```{raw} jupyter +
+ + QuantEcon + +
+``` + +# A Robust LQ Bewley Model + +```{contents} Contents +:depth: 2 +``` + +```{index} single: Robust Bewley Model +``` + +```{index} single: Observational Equivalence; heterogeneous agents +``` + +## Overview + +This lecture embeds the Bewley economy of {doc}`lq_bewley_complete_markets` in the robust permanent income framework of {doc}`lq_robust_smoothing`. + +It is the last of four lectures on the LQ permanent income model. + +The result is a family of economies in which consumers disagree about the model generating their income, yet behave identically. + +Using the observational-equivalence theorem of {cite:t}`HST_1999`, we show + +- how a continuum of consumers $i$ can differ in their robustness parameters $\sigma_i \leq 0$ and their discount factors $\beta_i$, provided each pair $(\sigma_i,\beta_i)$ lies on an observational-equivalence locus +- how every such consumer chooses the **same consumption-saving rule** as a benchmark $(\sigma = 0, \beta)$ agent who fully trusts the endowment process +- how the equilibrium interest rate $R = \beta^{-1}$ and all aggregate and cross-section dynamics therefore coincide with those of the plain-vanilla Bewley model +- how, despite all of this, distinct $(\sigma_i,\beta_i)$ agents hold genuinely different subjective models of their non-financial income + +The economy is a pure endowment economy, so there is no physical capital and investment plays no role. + +We read {doc}`lq_robust_smoothing` as a prerequisite and carry over its notation. + +```{note} +As in {doc}`lq_robust_smoothing`, $w_{t+1}$ is the baseline shock and $v_{t+1}$ a distortion to its conditional mean, $\sigma \le 0$ is the robustness parameter, $\eta_1$ and $\eta_2$ are the standard deviations of the permanent and transitory endowment shocks, and $a_t$ denotes net assets. +``` + +Let's begin with some imports. + +```{code-cell} ipython3 +import numpy as np +import matplotlib.pyplot as plt +``` + +## Mapping the Bewley economy into the HST framework + +We specialise the robust model of {doc}`lq_robust_smoothing` to $\lambda = \delta_h = 0$, so there are no habits and no durable goods, and to a pure endowment economy with no physical capital, $k_t = 0$. + +In this case services equal consumption, $s_t = c_t$. + +The only traded security is a one-period risk-free bond, and $a_t$ denotes the household's net asset position, so that positive $a_t$ is wealth. + +The endowment process follows the state-space representation + +$$ +\begin{aligned} +z_{t+1} &= \check{A}\, z_t + \check{C}\, w_{t+1} \\ +y_t &= \check{G}\, z_t +\end{aligned} +$$ (eq:rbew-endowment) + +with the two-factor specification $y_t = z_{1t}+z_{2t}$, $\check A = \mathrm{diag}(1,0)$ and $\check C = \mathrm{diag}(\eta_1,\eta_2)$. + +The household's augmented state vector is $x_t = [a_t,\; z_t^\top]^\top$, and the law of motion of {doc}`lq_robust_smoothing` specialises to + +$$ +\begin{pmatrix} a_{t+1} \\ z_{t+1} \end{pmatrix} += +\underbrace{\begin{pmatrix} R & R\check{G} \\ 0 & \check{A} \end{pmatrix}}_{A} +\begin{pmatrix} a_t \\ z_t \end{pmatrix} ++ +\underbrace{\begin{pmatrix} -R \\ 0 \end{pmatrix}}_{B} +c_t ++ +\underbrace{\begin{pmatrix} 0 \\ \check{C} \end{pmatrix}}_{C} +(w_{t+1} + v_{t+1}) +$$ (eq:rbew-law) + +The objective is $\mathbb{E}_0 \sum_{t=0}^\infty \beta^t\bigl[-(c_t - b)^2/2\bigr]$, which is the HST criterion with $\sigma = 0$ and a constant bliss level $b_t \equiv b$. + +The robust Bellman equation with $\sigma = 0$ therefore reduces exactly to the LQ problem of {doc}`lq_permanent_income`, confirming that the HST framework nests the Bewley model. + +## The robustness scalar for this economy + +Everything about the endowment process that matters for robustness is summarised by the scalar $\alpha^2$ derived in {doc}`lq_robust_smoothing`. + +For the two-factor endowment it is + +$$ +\alpha^2 = \eta_1^2 + (1-\beta)^2\,\eta_2^2 +$$ (eq:rbew-alpha) + +This is the variance of the consumption innovation $h\,w_{t+1}$, where $h = (1-\beta)\check G(I-\beta\check A)^{-1}\check C = \begin{pmatrix}\eta_1 & (1-\beta)\eta_2\end{pmatrix}$. + +In the present setting $\alpha^2$ has a second, equally concrete meaning that we met in {doc}`lq_bewley_complete_markets`. + +Because individual consumption is a random walk with innovation variance $\alpha^2$, the cross-section variance of consumption among agents of age $t$ who started from a common initial condition is exactly $t\,\alpha^2$. + +So the same scalar that governs how much a consumer's robustness concern bites also governs how fast the Bewley cross-section fans out. + +## The Bewley observational-equivalence locus + +Applying the observational-equivalence theorem {prf:ref}`thm-rcs-oe1` of {doc}`lq_robust_smoothing` at the equilibrium interest rate $R = \beta^{-1}$ gives the **Bewley observational-equivalence locus** + +$$ +\hat\beta(\sigma) = \beta + \frac{\sigma\,\alpha^2\,\beta}{1-\beta} +$$ (eq:rbew-locus) + +For $\sigma < 0$ we have $\hat\beta(\sigma) < \beta$. + +An agent with the pair $(\sigma, \hat\beta(\sigma))$ is more concerned about model misspecification, because $\sigma$ is lower, but also more impatient, because $\hat\beta$ is lower. + +The two forces cancel exactly, leaving the consumption decision rule unchanged. + +The locus is admissible only down to the breakdown point of {doc}`lq_robust_smoothing`, + +$$ +\underline\sigma = -\frac{(1-\beta)^2}{\alpha^2}, +\qquad\text{at which}\qquad +\hat\beta(\underline\sigma) = \beta^2 +$$ (eq:rbew-breakdown) + +Below $\underline\sigma$ the individual robust control problem has no solution, so there is no economy to describe. + +## Equilibrium with heterogeneous types + +We can now populate the economy with a continuum of types that differ in their concern for robustness. + +````{prf:proposition} A robust Bewley equilibrium +:label: prop-rbew-types + +Let each agent $i$ in the unit interval be indexed by a robustness parameter $\sigma_i \in (\underline\sigma, 0]$, distributed according to any distribution $\Phi$, and let agent $i$ have discount factor + +$$ +\beta_i = \hat\beta(\sigma_i) = \beta + \frac{\sigma_i\,\alpha^2\,\beta}{1-\beta} +$$ (eq:rbew-types) + +so that every pair $(\sigma_i,\beta_i)$ lies on the locus {eq}`eq:rbew-locus`. + +Then + +1. every agent's optimal consumption plan is identical to that of the plain-vanilla $(\sigma = 0,\, \beta)$ agent, +2. the equilibrium gross interest rate is $R = \beta^{-1}$, independently of $\Phi$, and +3. the aggregate and cross-section dynamics coincide with those of the benchmark Bewley economy of {doc}`lq_bewley_complete_markets`. +```` + +````{prf:proof} +By {prf:ref}`thm-rcs-oe1`, an agent with parameters $(\sigma_i, \hat\beta(\sigma_i))$ facing gross interest rate $R = \beta^{-1}$ chooses the same consumption-saving rule as the benchmark $(0,\beta)$ agent. + +This holds agent by agent and does not require the $\sigma_i$ to be equal, which gives part 1. + +Since all individual rules coincide with the benchmark rule, the goods-market clearing condition $\int c_t^i\, di = Y$ and the bond-market condition $\int a_t^i\, di = 0$ are the benchmark conditions, so they are satisfied at $R = \beta^{-1}$ for exactly the reason given in {doc}`lq_bewley_complete_markets`. + +Because market clearing never refers to $\Phi$, the equilibrium interest rate does not either, which gives part 2. + +Part 3 follows because aggregate and cross-section objects are integrals of individual paths, and the individual paths are the benchmark paths. +```` + +The distribution $\Phi$ of robustness types is therefore completely unidentified by quantity data. + +An econometrician who observes $\{c_t^i, a_t^i\}$ for every agent and every date cannot tell whether the economy is populated entirely by $\sigma_i = 0$ agents, entirely by $\sigma_i$ near $\underline\sigma$ agents, or by any mixture. + +## Where the agents genuinely differ + +Agents on the locus are indistinguishable in what they *do* but not in what they *believe*. + +An agent with $\sigma_i < 0$ applies a worst-case distortion $v_{t+1}^i = K(\sigma_i,\beta_i)\,\mu_{s,t}^i$ to its conditional expectations, while an agent with $\sigma_i = 0$ takes the approximating model at face value. + +From {doc}`lq_robust_smoothing`, the worst-case law for agent $i$'s marginal utility is + +$$ +\mu_{s,t+1}^i = \zeta_i\, \mu_{st}^i + \alpha\, w_{t+1}, +\qquad +\zeta_i = \frac{\beta}{\beta_i} = \left[1 + \frac{\sigma_i\alpha^2}{1-\beta}\right]^{-1} > 1 +$$ (eq:rbew-zeta) + +With $\lambda = \delta_h = 0$ and a constant bliss point we have $\mu_{st} = b - c_t$, so agent $i$'s **worst-case expected consumption path** is + +$$ +\hat{\mathbb{E}}_t\, c_{t+h}^i = b - \zeta_i^{\,h}\,(b - c_t) +$$ (eq:rbew-beliefs) + +Under the approximating model, by contrast, consumption is a martingale, $\mathbb{E}_t c_{t+h} = c_t$ for every $h$. + +Equation {eq}`eq:rbew-beliefs` says that a robust agent below its bliss point expects, under its worst-case model, that consumption will *drift away* from bliss at the geometric rate $\zeta_i$. + +The more robust the agent, the faster the drift it guards against, and the more precautionary saving it does. + +That extra saving is exactly offset by the lower $\beta_i$, which is why the realized path is the same for all types. + +## Computation + +We use the calibration of the preceding lectures. + +```{code-cell} ipython3 +β = 0.95 # benchmark discount factor, R = 1/β +η1 = 0.15 # std of permanent shock +η2 = 0.30 # std of transitory shock +b = 1.0 # bliss point + +R = 1 / β +h = np.array([η1, (1 - β) * η2]) # consumption innovation loadings +α2 = h @ h +α = np.sqrt(α2) +σ_lo = -(1 - β)**2 / α2 # breakdown point, eq:rbew-breakdown + +print(f"α^2 = {α2:.6f} α = {α:.6f}") +print(f"breakdown σ̲ = {σ_lo:.6f}, where β̂ = {β + σ_lo * α2 * β / (1 - β):.6f}" + f" (β² = {β**2:.6f})") +``` + +Next we build a set of types spread across the admissible range and record what distinguishes them. + +We reuse the detection error probability of {doc}`lq_robust_smoothing` to report how hard each type's worst-case model would be to detect in a sample of $T = 40$ quarters. + +```{code-cell} ipython3 +def worst_case_persistence(σ, β, α2): + "Worst-case persistence ζ(σ) of marginal utility, eq:rbew-zeta." + return 1 / (1 + σ * α2 / (1 - β)) + + +def simulate_paths(ζ, α, T, n_paths, seed): + "Simulate n_paths draws of μ_{t+1} = ζ μ_t + α w_{t+1} from μ_0 = 0." + rng = np.random.default_rng(seed) + paths = np.zeros((n_paths, T + 1)) + shocks = rng.standard_normal((n_paths, T)) + for t in range(T): + paths[:, t + 1] = ζ * paths[:, t] + α * shocks[:, t] + return paths + + +def log_likelihood_ratio(paths, ζ, α): + "Return log p_worst(path) - log p_approx(path)." + lag, lead = paths[:, :-1], paths[:, 1:] + return 0.5 * (np.sum(((lead - lag) / α)**2, axis=1) + - np.sum(((lead - ζ * lag) / α)**2, axis=1)) + + +def detection_error_probability(ζ, α, T=40, n_paths=10_000, seed=1234): + "Finite-sample DEP for the approximating and worst-case scalar laws." + if np.isclose(ζ, 1.0): + return 0.5 + approx = simulate_paths(1.0, α, T, n_paths, seed) + worst = simulate_paths(ζ, α, T, n_paths, seed + 1) + return 0.5 * (np.mean(log_likelihood_ratio(worst, ζ, α) < 0) + + np.mean(log_likelihood_ratio(approx, ζ, α) > 0)) +``` + +```{code-cell} ipython3 +σ_types = np.array([0.0, 0.3, 0.6, 0.9]) * σ_lo +β_types = β + σ_types * α2 * β / (1 - β) +ζ_types = worst_case_persistence(σ_types, β, α2) +dep_types = np.array([detection_error_probability(ζ, α) for ζ in ζ_types]) + +print(f"{'σ_i':>10}{'β_i':>10}{'ζ_i':>10}{'DEP':>8}") +for σ_i, β_i, ζ_i, dep in zip(σ_types, β_types, ζ_types, dep_types): + print(f"{σ_i:10.4f}{β_i:10.4f}{ζ_i:10.4f}{dep:8.3f}") +``` + +These four types differ substantially in patience and in the pessimism of their worst-case model. + +We now confirm that they nonetheless behave identically. + +Each type faces the same shocks and starts from the same initial consumption, and each consumes according to the benchmark rule $c_{t+1} = c_t + h\,w_{t+1}$. + +```{code-cell} ipython3 +T = 60 +rng = np.random.default_rng(42) +shocks = rng.standard_normal((T, 2)) # common shocks for all types + +c_paths = np.zeros((len(σ_types), T + 1)) +for i in range(len(σ_types)): + for t in range(T): + c_paths[i, t + 1] = c_paths[i, t] + h @ shocks[t] + +print("max absolute difference across types:" + f" {np.abs(c_paths - c_paths[0]).max():.2e}") +``` + +The paths agree exactly, which is {prf:ref}`prop-rbew-types` part 1 in action. + +The next figure contrasts what the types do with what they believe. + +```{code-cell} ipython3 +--- +mystnb: + figure: + caption: | + Same actions, different beliefs. Left: realized consumption paths for + four robustness types facing common shocks; the curves lie exactly on + top of one another. Right: each type's worst-case expected consumption + path from a common date, against the flat martingale forecast of the + approximating model. + name: fig-rbew-beliefs +--- +fig, axes = plt.subplots(1, 2, figsize=(11, 4)) + +for i, σ_i in enumerate(σ_types): + axes[0].plot(c_paths[i], lw=3 - 0.6 * i, alpha=0.9, + label=rf'$\sigma_i={σ_i:.3f}$') +axes[0].set_xlabel('$t$') +axes[0].set_ylabel('$c_t$') +axes[0].set_title('realized consumption') +axes[0].legend() + +horizons = np.arange(41) +c_now = c_paths[0, 20] +axes[1].axhline(c_now, color='k', linestyle=':', lw=1.2, + label='approximating model') +for σ_i, ζ_i in zip(σ_types, ζ_types): + if σ_i == 0.0: + continue + axes[1].plot(horizons, b - ζ_i**horizons * (b - c_now), lw=2, + label=rf'worst case, $\sigma_i={σ_i:.3f}$') +axes[1].set_xlabel('horizon $h$') +axes[1].set_ylabel(r'$\hat{\mathbb{E}}_t\,c_{t+h}$') +axes[1].set_title('expected consumption under each belief') +axes[1].legend() + +fig.tight_layout() +plt.show() +``` + +The left panel of {numref}`fig-rbew-beliefs` shows four curves drawn on top of one another. + +The right panel shows that the same four agents expect very different futures. + +The $\sigma_i = 0$ agent expects consumption to stay where it is. + +Every robust agent guards against a future in which consumption drifts away from the bliss point, and the drift is faster the more robust the agent. + +Finally we check part 3 of {prf:ref}`prop-rbew-types`, that the cross-section behaves as in the benchmark Bewley economy. + +We simulate a large population in which each agent draws its own type and its own shocks, and compare the cross-section variance of consumption to the benchmark prediction $t\,\alpha^2$. + +```{code-cell} ipython3 +n_agents, T = 20_000, 40 +rng = np.random.default_rng(1234) + +# each agent draws a robustness type; types do not affect behaviour +σ_i = rng.uniform(σ_lo, 0.0, size=n_agents) + +shocks = rng.standard_normal((n_agents, T, 2)) +c = np.zeros((n_agents, T + 1)) +c[:, 1:] = np.cumsum(shocks @ h, axis=1) + +print(f"{'t':>5}{'cross-section var':>20}{'t·α²':>12}") +for t in [10, 20, 30, 40]: + print(f"{t:5d}{c[:, t].var():20.5f}{t * α2:12.5f}") +``` + +The cross-section variance grows linearly at rate $\alpha^2$, exactly as in {doc}`lq_bewley_complete_markets`, and the distribution of robustness types leaves no trace in the data. + +## Concluding remarks + +We embedded the single-agent robust permanent income model in a Bewley equilibrium with a continuum of agents who differ in how much they distrust their income process. + +Provided each agent's pair $(\sigma_i,\beta_i)$ lies on the observational-equivalence locus {eq}`eq:rbew-locus`, every agent chooses the benchmark consumption rule. + +The equilibrium interest rate $R = \beta^{-1}$, the aggregate dynamics, and the linear growth of the cross-section variance of consumption are all inherited unchanged from the plain-vanilla Bewley model of {doc}`lq_bewley_complete_markets`. + +This is a strong non-identification result. + +Quantity data pin down the decision rule, but they cannot decompose it into a degree of impatience and a degree of concern about misspecification. + +What the agents do *not* share is their view of the world: each robust type acts as if its consumption were about to drift away from bliss at its own rate $\zeta_i$. + +Two routes lead out of this observational equivalence. + +One is to look at asset prices, which do distinguish the $(\sigma,\hat\beta)$ pairs, and which are studied in {doc}`robust_permanent_income`. + +The other is to bound the plausible range of $\sigma$ statistically, using the detection error probabilities of {doc}`lq_robust_smoothing`. + +## Exercises + +```{exercise-start} +:label: rbew_ex1 +``` + +This exercise asks you to carry out the translation from the benchmark Bewley economy into HST notation. + +Specialise the robust-control setup to the no-habit, no-capital LQ Bewley environment, so $\lambda = \delta_h = 0$ and $k_t = 0$, and let the endowment follow the two-factor model. + +1. Write the household state as $x_t = [a_t, z_t^\top]^\top$, where $a_t$ is net assets, and derive the matrices $(A, B, C)$ in {eq}`eq:rbew-law`. + +2. Show that when $\sigma = 0$ the Bellman problem coincides with the LQ permanent-income problem of {doc}`lq_permanent_income`. + +3. HST define $\alpha^2 = \nu^\top\nu$ with $\nu^\top = M_s C$, where $\mu_{st} = M_s x_t$. Compute $M_s$ for this economy and verify that this route delivers the same $\alpha^2$ as {eq}`eq:rbew-alpha`. + +```{exercise-end} +``` + +```{solution-start} rbew_ex1 +:class: dropdown +``` + +Here is one solution. + +1. With budget law $a_{t+1} = R(a_t + y_t - c_t)$, $y_t = \check G z_t$ and $z_{t+1} = \check A z_t + \check C w_{t+1}$, stacking gives + +$$ +\begin{pmatrix} a_{t+1} \\ z_{t+1} \end{pmatrix} += +\underbrace{\begin{pmatrix} R & R\check G \\ 0 & \check A \end{pmatrix}}_{A} +\begin{pmatrix} a_t \\ z_t \end{pmatrix} ++ +\underbrace{\begin{pmatrix} -R \\ 0 \end{pmatrix}}_{B} c_t ++ +\underbrace{\begin{pmatrix} 0 \\ \check C \end{pmatrix}}_{C} w_{t+1} . +$$ + + The sign of $B$ is negative because higher $c_t$ reduces asset accumulation. + +2. At $\sigma = 0$ the minimizing agent is absent, the distortion term drops out of the Bellman equation, and the objective is $\mathbb{E}_0\sum \beta^t[-(c_t-b)^2/2]$ subject to a linear law of motion. + + That is precisely the LQ permanent-income problem. + +3. With $\lambda = \delta_h = 0$ and constant bliss, $\mu_{st} = b - c_t$, and the optimal rule is + +$$ +c_t = (1-\beta)\bigl[a_t + \check G(I-\beta\check A)^{-1} z_t\bigr] + \text{constant}, +$$ + + so $M_s = -(1-\beta)\begin{pmatrix}1 & \check G(I-\beta\check A)^{-1}\end{pmatrix}$ up to the constant. + + Since $C = [0;\ \check C]$, the first column of $M_s$ is annihilated and + +$$ +\nu^\top = M_s C = -(1-\beta)\check G(I-\beta\check A)^{-1}\check C = -h . +$$ + + Hence $\alpha^2 = \nu^\top\nu = hh^\top = \eta_1^2 + (1-\beta)^2\eta_2^2$, matching {eq}`eq:rbew-alpha`. + + The minus sign is immaterial because only $\alpha^2$ appears. + +```{solution-end} +``` + +```{exercise-start} +:label: rbew_ex2 +``` + +This exercise works through the equilibrium logic of {prf:ref}`prop-rbew-types`. + +Fix a benchmark pair $(\beta, \sigma = 0)$ with $R = \beta^{-1}$ and let a unit interval of consumers be indexed by $i$, with type $\sigma_i \in (\underline\sigma, 0]$ and discount factor $\beta_i = \hat\beta(\sigma_i)$ from {eq}`eq:rbew-locus`. + +1. Use {prf:ref}`thm-rcs-oe1` of {doc}`lq_robust_smoothing` to show that each type has the same consumption rule as the benchmark $(\beta, 0)$ agent. + +2. Show that goods- and bond-market clearing imply the same equilibrium interest rate $R = \beta^{-1}$ as in the plain-vanilla Bewley model, whatever the distribution of types. + +3. Explain why agents can be observationally equivalent in quantities while holding different worst-case subjective models. + +4. Suppose instead that agents share a common $\beta$ but differ in $\sigma_i$, so that they are *not* on the locus. Explain why $R = \beta^{-1}$ would then generally fail to clear the bond market. + +```{exercise-end} +``` + +```{solution-start} rbew_ex2 +:class: dropdown +``` + +Here is one solution. + +1. {prf:ref}`thm-rcs-oe1` says that if $(\sigma_i,\beta_i)$ satisfies $\beta_i = \beta + \sigma_i\alpha^2\beta/(1-\beta)$, then type $i$ chooses the same decision rule as the benchmark agent, so all types share the policy function $c_t = \mathcal{C}(a_t, z_t)$. + +2. Since all individual rules coincide with the benchmark rule, aggregating over $i$ reproduces the benchmark market-clearing conditions, which hold at $R = \beta^{-1}$. + + The type distribution never enters, so it cannot affect the equilibrium rate. + +3. Observational equivalence is a statement about quantities generated by optimal rules. + + The minimizing feedback $K(\sigma_i,\beta_i)$ still differs across types, so the agents attach different worst-case conditional means to the same shock process while making identical choices. + +4. Off the locus the impatience offset is missing. + + An agent with $\sigma_i < 0$ and discount factor $\beta$ has a precautionary motive that is not cancelled, so it wants to save more than the benchmark agent at $R = \beta^{-1}$. + + With positive net demand for bonds in the aggregate, the equilibrium interest rate must fall below $\beta^{-1}$ to clear the market, and the equilibrium then depends on the whole distribution of types. + +```{solution-end} +``` + +```{exercise-start} +:label: rbew_ex3 +``` + +This exercise separates quantities from beliefs for two individual agents. + +Consider agents $a$ and $b$ with $\sigma^a < \sigma^b \leq 0$, both on the locus {eq}`eq:rbew-locus`. + +1. Show that the two agents have the same consumption innovation $h\,w_{t+1}$. + +2. Show that if they start from the same $(a_t, z_t)$ and observe the same shock $w_{t+1}$, their next-period consumption and assets coincide. + +3. Using {eq}`eq:rbew-beliefs`, compute the ratio of the two agents' worst-case forecasts of $b - c_{t+h}$ and show that it grows geometrically in $h$. + +4. Summarise what is and is not identified by data on quantities alone. + +```{exercise-end} +``` + +```{solution-start} rbew_ex3 +:class: dropdown +``` + +Here is one solution. + +1. Both pairs lie on {eq}`eq:rbew-locus`, so by {prf:ref}`thm-rcs-oe1` both use the benchmark rule and hence the same innovation vector $h$. + +2. With a common state and a common shock, both agents apply the same policy function and the same law of motion, so $c_{t+1}^a = c_{t+1}^b$ and $a_{t+1}^a = a_{t+1}^b$. + +3. From {eq}`eq:rbew-beliefs`, $\hat{\mathbb{E}}_t (b - c_{t+h}^j) = \zeta_j^{\,h}(b-c_t)$, so the ratio is $(\zeta_a/\zeta_b)^h$. + + Since $\sigma^a < \sigma^b$ implies $\zeta_a > \zeta_b$, the ratio grows geometrically: the two agents' beliefs diverge without bound as the horizon lengthens, even though their actions never differ at all. + +4. Quantities identify the equilibrium decision rule, and hence the single combination of parameters that appears in it. + + They do not identify the decomposition of that rule into impatience $\beta_i$ and robustness $\sigma_i$ along the locus. + +```{solution-end} +``` + +```{exercise-start} +:label: rbew_ex4 +``` + +This exercise asks how much belief heterogeneity is statistically plausible. + +Restrict attention to types whose worst-case model has a detection error probability of at least $0.2$ in a sample of $T = 40$. + +1. Find the most robust admissible type $\sigma^{\min}$ by bisection. + +2. For that type, report $\beta_i$, $\zeta_i$, and the horizon $h$ at which its worst-case forecast of $b - c_{t+h}$ is twice the approximating model's forecast. + +3. Repeat part 1 with $T = 160$ and comment on what a longer sample does to the plausible amount of belief heterogeneity. + +```{exercise-end} +``` + +```{solution-start} rbew_ex4 +:class: dropdown +``` + +Here is one solution. + +The approximating model forecasts $b - c_{t+h} = b - c_t$ for every $h$, while the worst-case forecast is $\zeta_i^h(b-c_t)$, so the doubling horizon solves $\zeta_i^h = 2$. + +```{code-cell} ipython3 +def σ_for_target_dep(target, T, β, α2, tol=1e-5): + """ + Find σ ∈ (σ̲, 0) with DEP(σ) = target by bisection. + + Returns None if the DEP never falls to the target on the admissible + range, in which case the breakdown point is the binding constraint. + """ + α_loc = np.sqrt(α2) + lo, hi = 0.999 * (-(1 - β)**2 / α2), 0.0 + + def dep_at(σ): + return detection_error_probability( + worst_case_persistence(σ, β, α2), α_loc, T=T) + + if dep_at(lo) > target: + return None + + while hi - lo > tol: + mid = 0.5 * (lo + hi) + if dep_at(mid) < target: + lo = mid + else: + hi = mid + return 0.5 * (lo + hi) + + +for T in [40, 160]: + σ_min = σ_for_target_dep(0.2, T, β, α2) + if σ_min is None: + σ_min = 0.999 * σ_lo # breakdown binds before detectability + note = ' (breakdown point binds)' + else: + note = '' + ζ_min = worst_case_persistence(σ_min, β, α2) + print(f"T = {T:>3}: σ_min = {σ_min:.5f} β_i = {β / ζ_min:.5f} " + f"ζ_i = {ζ_min:.5f} doubling horizon = " + f"{np.log(2) / np.log(ζ_min):.1f} quarters{note}") +``` + +At $T = 40$ the DEP never falls to $0.2$ within the admissible range, so the most robust plausible type is the one at the breakdown point itself. + +At $T = 160$ statistical detectability binds first and the plausible set of types shrinks sharply toward $\sigma = 0$. + +The doubling horizon lengthens correspondingly: with more data, only agents whose pessimism accumulates slowly remain statistically credible. + +```{solution-end} +``` + +## Related lectures + +- {doc}`lq_permanent_income` develops the standard LQ permanent income model. +- {doc}`lq_bewley_complete_markets` builds the benchmark Bewley economy whose equilibrium is reproduced here. +- {doc}`lq_robust_smoothing` derives the observational-equivalence theorems, the breakdown point, and the detection error probabilities used in this lecture. +- {doc}`robust_permanent_income` shows how asset prices break the observational equivalence. diff --git a/lectures/lq_robust_smoothing.md b/lectures/lq_robust_smoothing.md index b49efc82c..07c994fcb 100644 --- a/lectures/lq_robust_smoothing.md +++ b/lectures/lq_robust_smoothing.md @@ -36,50 +36,49 @@ kernelspec: This lecture studies a robust version of the LQ permanent income model due to {cite:t}`HST_1999` and {cite:t}`HansenSargent2008`. -It is the third of three lectures on the LQ permanent income model. +A consumer who distrusts his specification of the labor income process engages in a form of precautionary savings. + +This is the third of four lectures on the LQ permanent income model. It builds on {doc}`lq_permanent_income`, which develops the standard model, and {doc}`lq_bewley_complete_markets`, which studies its cross-section and market-structure implications. -A consumer who distrusts his specification of the labor income process engages in a form of precautionary savings. +The sequel, {doc}`lq_robust_bewley`, uses the results developed here to build a Bewley economy populated by consumers who differ in how much they distrust their income model. Our description of the model with concerns about robustness includes -- how (for quantities) a concern for robustness is observationally equivalent to an increase in - impatience +- how, for quantities, a concern for robustness is observationally equivalent to an increase in impatience - how the worst-case model that the consumer uses to shape his decision rule distorts the baseline model's endowment process toward greater persistence +- a **breakdown point** beyond which the robust control problem ceases to have a solution - a frequency-domain representation of the effects of concerns about misspecification of the endowment process - a detection-error-probability characterization of the amount of model uncertainty -The lecture concludes by combining the Bewley economy of {doc}`lq_bewley_complete_markets` with the robustness machinery. - -Using tools from {cite:t}`HansenSargent2008`, we show: - -- how a continuum of consumers $i$ who use identical decision rules can nevertheless differ in their robustness parameters $\sigma_i \leq 0$ and - their discount factors $\beta_i$, provided that the pair $(\sigma_i, \beta_i)$ lies on an observational-equivalence locus - derived below -- how every such consumer chooses the **same consumption-saving rule** as a baseline - plain-vanilla $(\sigma = 0, \beta)$ agent with no concerns about misspecification of the endowment process -- how the equilibrium interest rate $R = \beta^{-1}$ and all aggregate dynamics therefore - coincide with those of a benchmark Bewley model -- how distinct $(\sigma_i, \beta_i)$ agents act as if they have different subjective models of their non-financial income process - -We first present the HST model in its general form, which includes physical capital and investment $i_t$. - -When we return to the Bewley economy of {doc}`lq_bewley_complete_markets`, we specialise to a pure endowment economy with no capital, so investment plays no role there. +A recurring theme is that a single scalar $\alpha^2$, the variance of the innovation to the consumer's marginal utility, summarises everything about the endowment process that matters for robustness. Let's begin with some imports. ```{code-cell} ipython3 import numpy as np import matplotlib.pyplot as plt -from scipy.stats import norm - ``` ## A brief review We recall the essentials from {doc}`lq_permanent_income` and {doc}`lq_bewley_complete_markets`. +### Notation + +Because a robust decision maker guards against distortions to the mean of a shock, we need separate symbols for the shock and for the distortion. + +We therefore adopt the following conventions, which differ in three places from the two preceding lectures. + +```{note} +- $w_{t+1}$ is the baseline IID shock, as in {doc}`lq_permanent_income`, and $v_{t+1}$ is a **distortion** to its conditional mean, as in {doc}`robust_permanent_income`. +- $\sigma \leq 0$ is the **robustness parameter**. The standard deviations of the two endowment shocks, written $\sigma_1$ and $\sigma_2$ in the preceding lectures, are renamed $\eta_1$ and $\eta_2$ here so that $\sigma$ is free. +- $a_t$ denotes the consumer's **net assets**, equal to minus the debt $b_t$ of {doc}`lq_permanent_income`. This frees $b_t$ for the preference shifter of {cite:t}`HST_1999`. +``` + +### The model + A consumer with quadratic utility and discount factor $\beta$ faces the endowment process $$ @@ -87,18 +86,18 @@ $$ z_{t+1} &= \check{A}\, z_t + \check{C}\, w_{t+1} \\ y_t &= \check{G}\, z_t \end{aligned} -$$ (eq:rs-endowment) +$$ (eq:rcs-endowment) The optimal decision rule has a state-space representation in which the state is current consumption $c_t$ and the exogenous endowment state $z_t$: $$ \begin{aligned} c_{t+1} &= c_t + (1-\beta)\,\check{G}(I-\beta\check{A})^{-1}\check{C}\, w_{t+1} \\ -b_t &= \check{G}(I-\beta\check{A})^{-1} z_t - \frac{1}{1-\beta}\,c_t \\ +a_t &= \frac{1}{1-\beta}\,c_t - \check{G}(I-\beta\check{A})^{-1} z_t \\ y_t &= \check{G}\, z_t \\ z_{t+1} &= \check{A}\, z_t + \check{C}\, w_{t+1} \end{aligned} -$$ (eq:rs-crep) +$$ (eq:rcs-crep) We again use the two-factor endowment $y_t = z_{1t} + z_{2t}$, @@ -108,48 +107,88 @@ $$ \begin{pmatrix}1 & 0\\0 & 0\end{pmatrix} \begin{pmatrix}z_{1t}\\z_{2t}\end{pmatrix} + -\begin{pmatrix}\sigma_1 & 0\\0 & \sigma_2\end{pmatrix} +\begin{pmatrix}\eta_1 & 0\\0 & \eta_2\end{pmatrix} \begin{pmatrix}w_{1,t+1}\\w_{2,t+1}\end{pmatrix} -$$ (eq:pi-twofactor) +$$ (eq:rcs-twofactor) with $z_{1t}$ a permanent component and $z_{2t}$ a purely transitory component. -The following cell fixes the calibration used below. +### The consumption innovation -```{code-cell} ipython3 -# Parameters (as in the preceding lectures) -β = 0.95 # discount factor -σ1 = 0.15 # std of permanent shock -σ2 = 0.30 # std of transitory shock -``` +One scalar built from {eq}`eq:rcs-crep` will do all of the work below. -## A robust permanent income model +The first line of {eq}`eq:rcs-crep` says that consumption is a random walk whose innovation is $h\, w_{t+1}$, where -```{index} single: Robust Control; permanent income +$$ +h = (1-\beta)\,\check{G}(I-\beta\check{A})^{-1}\check{C} +$$ (eq:rcs-h) + +Define $\alpha^2$ to be the variance of that innovation, + +$$ +\alpha^2 = h h^\top += (1-\beta)^2\,\check{G}(I-\beta\check{A})^{-1}\check{C}\check{C}^\top(I-\beta\check{A}^\top)^{-1}\check{G}^\top +$$ (eq:rcs-alpha) + +For the two-factor endowment {eq}`eq:rcs-twofactor` we have $\check A = \mathrm{diag}(1,0)$, $\check C = \mathrm{diag}(\eta_1,\eta_2)$ and $\check G = \begin{pmatrix}1 & 1\end{pmatrix}$, so that $(I-\beta\check A)^{-1} = \mathrm{diag}\bigl((1-\beta)^{-1},1\bigr)$ and + +$$ +h = \begin{pmatrix}\eta_1 & (1-\beta)\eta_2\end{pmatrix}, +\qquad +\alpha^2 = \eta_1^2 + (1-\beta)^2\,\eta_2^2 +$$ (eq:rcs-alpha2) + +The permanent shock variance $\eta_1^2$ enters with coefficient $1$ because a unit permanent shock is *fully* capitalised into consumption. + +The transitory shock variance $\eta_2^2$ enters with the small coefficient $(1-\beta)^2$ because only its annuity value is consumed. + +This scalar does triple duty across the three lectures of this suite. + +```{note} +$\alpha^2$ is simultaneously + +- the variance of the consumption innovation in {doc}`lq_permanent_income`, +- the rate at which the cross-section variance of consumption grows with age in {doc}`lq_bewley_complete_markets`, and +- the quantity that, multiplied by $\sigma$, governs every robustness result in this lecture. + +{doc}`robust_permanent_income` writes the same object as $\theta^2$. ``` -```{index} single: Precautionary Savings; robustness +The following cell fixes the calibration used below. + +```{code-cell} ipython3 +β = 0.95 # discount factor, so R = 1/β +η1 = 0.15 # std of permanent shock +η2 = 0.30 # std of transitory shock + +R = 1 / β +α2 = η1**2 + (1 - β)**2 * η2**2 +α = np.sqrt(α2) + +print(f"α^2 = {α2:.6f}") +print(f" permanent η1^2 = {η1**2:.6f} " + f"({100 * η1**2 / α2:5.1f}% of α^2)") +print(f" transitory (1-β)^2 η2^2 = {(1 - β)**2 * η2**2:.6f} " + f"({100 * (1 - β)**2 * η2**2 / α2:5.1f}% of α^2)") ``` +Permanent shocks account for almost all of $\alpha^2$ in this calibration. + +## A robust permanent income model + ### Robustness and precautionary savings We now study a consumer who *distrusts* his specification of the stochastic process governing his labor income. -The model is due to {cite:t}`HST_1999` (HST), who estimated -it on US quarterly consumption and investment data. +The model is due to {cite:t}`HST_1999` (HST), who estimated it on US quarterly consumption and investment data. For a fuller treatment of the HST model and its asset-pricing implications, see {doc}`robust_permanent_income`. -A consumer who fears model misspecification engages in a form of **precautionary savings** that is -distinct from the usual precautionary motive (which requires a convex marginal utility). +A consumer who fears model misspecification engages in a form of **precautionary savings** that is distinct from the usual precautionary motive, which requires a convex marginal utility. -Here, the -precautionary motive arises because the consumer wants to protect against misspecification of the -**conditional means** of income shocks, and it operates even with quadratic preferences. +Here, the precautionary motive arises because the consumer wants to protect against misspecification of the **conditional means** of income shocks, and it operates even with quadratic preferences. -HST showed an important **observational equivalence** result: for quantities $(c_t, i_t)$ alone, -a concern for robustness is indistinguishable from an increase in impatience (a decrease in -$\beta$). +HST showed an important **observational equivalence** result: for quantities $(c_t, i_t)$ alone, a concern for robustness is indistinguishable from an increase in impatience, that is, a decrease in $\beta$. We develop this result carefully below. @@ -161,17 +200,16 @@ We develop this result carefully below. ```{index} single: Hansen Sargent Tallarini; model ``` -HST's model features a planner with preferences over consumption streams $\{c_t\}$, mediated -through **service streams** $\{s_t\}$. +HST's model features a planner with preferences over consumption streams $\{c_t\}$, mediated through **service streams** $\{s_t\}$. -Let $b$ be a preference shifter (utility bliss point). +Let $b$ be a preference shifter, or utility bliss point. The **Bellman equation for the robust planner** is $$ -x^\top P x - p = -\sup_c \inf_w \Bigl\{-(s-b)^2 + \beta\bigl(\theta (w^*)^\top w^* - \mathbb{E}\,(x^*)^\top P x^* - p\bigr)\Bigr\} -$$ (eq:income1) +\sup_c \inf_{v^*} \Bigl\{-(s-b)^2 + \beta\bigl(\theta\, (v^*)^\top v^* - \mathbb{E}\,(x^*)^\top P x^* - p\bigr)\Bigr\} +$$ (eq:rcs-bellman) subject to the household technology, capital accumulation, endowment dynamics, and the state law: @@ -182,70 +220,55 @@ h^* &= \delta_h h + (1-\delta_h) c \\ k^* &= \delta_k k + i \\ c + i &= \gamma k + d \\ \begin{pmatrix}d\\b\end{pmatrix} &= U z \\ -z^* &= A_{22} z + C_2(\epsilon^* + w^*) +z^* &= A_{22} z + C_2(w^* + v^*) \end{aligned} -$$ (eq:income1a) +$$ (eq:rcs-tech) -Here $^*$ denotes the next-period value; $c$ is consumption; $s$ is the scalar service measure; -$h$ is a habit stock; $k$ is the capital stock; $i$ is investment; $d$ is an endowment/technology -shock; $b$ is a **preference shock** (bliss-point shifter, distinct from the bond/debt variable -$b_t$ used above); $\epsilon^* \sim N(0,I)$ is the baseline shock; and -$w^*$ is a **distortion** to the conditional mean of $\epsilon^*$ chosen by a minimizing agent. +Here $^*$ denotes the next-period value; $c$ is consumption; $s$ is the scalar service measure; $h$ is a habit stock; $k$ is the capital stock; $i$ is investment; $d$ is an endowment shock; $b$ is a **preference shock**; $\gamma$ is the marginal product of capital; $w^* \sim N(0,I)$ is the baseline shock; and $v^*$ is a **distortion** to the conditional mean of $w^*$ chosen by a minimizing agent. -The penalty parameter $\theta > 0$ governs the consumer's concern about robustness. +The penalty parameter $\theta$ governs the consumer's concern about robustness. + +A large $\theta$ makes distortions expensive and so restrains the minimizing agent. We use the transformation $$ -\sigma = -\theta^{-1} \leq 0 -$$ +\sigma = -\theta^{-1} \leq 0, +\qquad \theta \in (0,\infty] +$$ (eq:rcs-sigma) -so $\sigma = 0$ corresponds to no robustness concern and $\sigma < 0$ to an increasing concern. +so that $\sigma = 0$, equivalently $\theta = \infty$, corresponds to no robustness concern, and $\sigma < 0$ to an increasing concern. -When $\lambda > 0$ and $\delta_h \in (0,1)$, the technology -{eq}`eq:income1a` accommodates **habit persistence** (positive $\lambda$) or durability. +When $\lambda > 0$ and $\delta_h \in (0,1)$, the technology {eq}`eq:rcs-tech` accommodates **habit persistence** or durability, and the stock $h_t$ is a geometric weighted average of current and past consumption. -The stock -$h_t$ is a geometric weighted average of current and past consumption. - -Equation $c_t + k_t = Rk_{t-1} + d_t$ with -$R = \delta_k + \gamma$ combines capital accumulation with a linear production technology. - -$R$ is -the physical gross return on capital. +Equation $c_t + k_t = R k_{t-1} + d_t$ with $R = \delta_k + \gamma$ combines capital accumulation with a linear production technology, so $R$ is the physical gross return on capital. Let $x_t^\top = [h_{t-1},\, k_{t-1},\, z_t^\top]$. -The state transition equations are: +The state transition equation is $$ -x_{t+1} = A\, x_t + B\, u_t + C(\epsilon_{t+1} + w_{t+1}) -$$ (eq:law0) +x_{t+1} = A\, x_t + B\, u_t + C(w_{t+1} + v_{t+1}) +$$ (eq:rcs-law) -where $u_t = c_t$ and $w_{t+1}$ is the distortion to the conditional mean of $\epsilon_{t+1}$. +where $u_t = c_t$ and $v_{t+1}$ is the distortion to the conditional mean of $w_{t+1}$. -HST estimated the model on U.S. quarterly data (1970Q1-1996Q3) using -nondurables plus services for consumption and durable consumption plus gross private investment for -investment. +HST estimated the model on US quarterly data from 1970Q1 to 1996Q3, using nondurables plus services for consumption and durable consumption plus gross private investment for investment. -Key estimates are summarised in the following table (reported in Appendix A of HST): +They imposed $\beta R = 1$ and $\delta_k = 0.975$, so $\gamma$ is pinned down once $\beta$ is estimated. -| Parameter | Habit | No Habit | +Two of their preference estimates are worth recording. + +| Parameter | Habit | No habit | |-----------|-------|----------| -| Risk-free rate | 0.025 | 0.025 | | $\beta$ | 0.997 | 0.997 | | $\delta_h$ | 0.682 | — | | $\lambda$ | 2.443 | 0 | -| $\alpha_1$ | 0.813 | 0.900 | -| $\alpha_2$ | 0.189 | 0.241 | -| $\phi_1$ | 0.998 | 0.995 | -| $\phi_2$ | 0.704 | 0.450 | | $2 \times \log L$ | 779.05 | 762.55 | -HST imposed $\beta R = 1$ and $\delta_k = 0.975$, so $\gamma$ is pinned down once $\beta$ is -estimated. +At a quarterly frequency, $\beta = 0.997$ implies an annual real interest rate of $\beta^{-4} - 1 \approx 1.2\%$. -An annual real interest rate of 2.5% corresponds to $\beta = 0.997$. +The remaining estimated parameters govern the exogenous $d_t$ and $b_t$ processes and are reported in Appendix A of {cite:t}`HST_1999`. ### Solution when $\sigma = 0$ @@ -253,9 +276,9 @@ When $\sigma = 0$ the objective reduces to $$ \mathbb{E}_0\sum_{t=0}^{\infty}\beta^t\bigl\{-(s_t - b_t)^2\bigr\} -$$ (eq:income5) +$$ (eq:rcs-obj) -Formulating a Lagrangian and deriving first-order conditions yields: +Forming a Lagrangian and deriving first-order conditions yields $$ \begin{aligned} @@ -264,18 +287,15 @@ $$ \mu_{ht} &= \beta \mathbb{E}_t[\delta_h \mu_{h,t+1} - \lambda \mu_{s,t+1}] \\ \mu_{ct} &= \beta R\, \mathbb{E}_t\mu_{c,t+1} \end{aligned} -$$ (eq:foc) +$$ (eq:rcs-foc) -Here $\mu_{st}$ is the **marginal valuation of consumption services**, which summarises the -endogenous state variables $h_{t-1}$ and $k_{t-1}$. +Here $\mu_{st}$ is the **marginal valuation of consumption services**, which summarises the endogenous state variables $h_{t-1}$ and $k_{t-1}$. -Equation {eq}`eq:foc` (last line) implies -$\mathbb{E}_t\mu_{c,t+1} = (\beta R)^{-1}\mu_{ct}$, so $\mu_{st}$ is a martingale -when $\beta R = 1$: +The last line of {eq}`eq:rcs-foc` implies $\mathbb{E}_t\mu_{c,t+1} = (\beta R)^{-1}\mu_{ct}$, so $\mu_{st}$ is a martingale when $\beta R = 1$: $$ -\mu_{st} = \mu_{s,t-1} + \nu^\top \epsilon_t -$$ (eq:martingale) +\mu_{st} = \mu_{s,t-1} + \nu^\top w_t +$$ (eq:rcs-martingale) for some vector $\nu$. @@ -284,108 +304,111 @@ Solving forward and substituting gives $$ \mu_{st} = \Psi_1 k_{t-1} + \Psi_2 h_{t-1} + \Psi_3\sum_{j=0}^{\infty} R^{-j} \mathbb{E}_t b_{t+j} + \Psi_4\sum_{j=0}^{\infty} R^{-j} \mathbb{E}_t d_{t+j} -$$ (eq:income10) +$$ (eq:rcs-mus) where $$ \Psi_1 = -(1+\lambda)R(1-R^{-2}\beta^{-1})\!\left[\frac{1-R^{-1}\tilde\delta_h}{1-R^{-1}\tilde\delta_h+\lambda(1-\tilde\delta_h)}\right], \quad \Psi_4 = R^{-1}\Psi_1 -$$ (eq:income100a) +$$ (eq:rcs-psi) and $\tilde\delta_h = (\delta_h + \lambda)/(1+\lambda)$. -In the widely-studied special case $\lambda = \delta_h = 0$, so $s_t = c_t$ and -$\mu_{st} = b_t - c_t$, the marginal propensity to consume out of **non-human wealth** $Rk_{t-1}$ -equals that out of **human wealth** $\sum_{j=0}^{\infty}R^{-j}\mathbb{E}_t d_{t+j}$, a well-known feature of -the LQ model. +In the widely-studied special case $\lambda = \delta_h = 0$, we have $s_t = c_t$ and $\mu_{st} = b_t - c_t$, and the marginal propensity to consume out of **non-human wealth** $Rk_{t-1}$ equals that out of **human wealth** $\sum_{j=0}^{\infty}R^{-j}\mathbb{E}_t d_{t+j}$, a well-known feature of the LQ model. -The formula for $\mu_{st}$ can be written as $\mu_{st} = M_s x_t$ where $x_t$ follows {eq}`eq:law0`. +The formula for $\mu_{st}$ can be written as $\mu_{st} = M_s x_t$ where $x_t$ follows {eq}`eq:rcs-law`. It follows that $$ \nu^\top = M_s C, \qquad \alpha = \sqrt{\nu^\top \nu} = \sqrt{M_s C C^\top M_s^\top} -$$ (eq:hsoffset2) +$$ (eq:rcs-nu) + +This $\alpha$ is the same scalar we met in {eq}`eq:rcs-alpha`. -The scalar $\alpha$ plays a central role in the observational equivalence result below. +To see why, set $\lambda = \delta_h = 0$ and hold $b_t$ fixed, so that $\mu_{st} = b - c_t$ and the innovation to $\mu_{st}$ is minus the innovation to $c_t$. -### Observational equivalence +Hence $\nu^\top = -h$ and $\alpha^2 = \nu^\top\nu = h h^\top$, exactly as in {eq}`eq:rcs-alpha`. + +The sign is irrelevant because only $\alpha^2$ ever appears. + +## Observational equivalence ```{index} single: Observational Equivalence; Theorem 1 ``` HST state an observational-equivalence theorem. -````{prf:theorem} Observational Equivalence, I -:label: thm-lqcs-oe1 +````{prf:theorem} Observational equivalence, I +:label: thm-rcs-oe1 Fix all parameters except $(\sigma, \beta)$ and suppose $\beta R = 1$ when $\sigma = 0$. -There exists $\underline\sigma < 0$ such that for any -$\sigma \in (\underline\sigma, 0)$, the optimal consumption-investment plan for $(0,\beta)$ is also -chosen by a robust decision maker with parameters $(\sigma, \hat\beta(\sigma))$, where +There exists $\underline\sigma < 0$ such that for any $\sigma \in (\underline\sigma, 0)$, the optimal consumption-investment plan for $(0,\beta)$ is also chosen by a robust decision maker with parameters $(\sigma, \hat\beta(\sigma))$, where $$ \hat\beta(\sigma) = \frac{1}{R} + \frac{\sigma\alpha^2}{R-1} -$$ (eq:obseq) += \beta + \frac{\sigma\alpha^2\beta}{1-\beta} +$$ (eq:rcs-oe) and $\hat\beta(\sigma) < \beta$. ```` -Since $R > 1$ and $\alpha^2 > 0$, a more negative $\sigma$ (stronger robustness -concern) lowers $\hat\beta$. +The second equality in {eq}`eq:rcs-oe` uses $R = \beta^{-1}$ and will be the form we use in computations. + +Since $R > 1$ and $\alpha^2 > 0$, a more negative $\sigma$, meaning a stronger robustness concern, lowers $\hat\beta$. A robust consumer wants to save more because his alter ego, a utility-minimizing agent, makes future income look worse than the approximating model predicts. A lower discount factor makes a consumer less patient and therefore reduces saving. -When these two forces are balanced according to {eq}`eq:obseq`, consumption plans are identical across $(\sigma, \hat\beta(\sigma))$ pairs. +When these two forces are balanced according to {eq}`eq:rcs-oe`, consumption plans are identical across $(\sigma, \hat\beta(\sigma))$ pairs. ````{prf:proof} When $\beta R = 1$ and $\sigma = 0$, the marginal utility $\mu_{st}$ obeys the martingale $$ -\mu_{st} = \mu_{s,t-1} + \alpha\,\tilde\epsilon_t -$$ (eq:reversee1) +\mu_{st} = \mu_{s,t-1} + \alpha\,\tilde w_t +$$ (eq:rcs-scalar-approx) -where $\tilde\epsilon_t$ is scalar IID with mean zero and unit variance. +where $\tilde w_t$ is scalar IID with mean zero and unit variance. -Activating a concern about robustness ($\sigma < 0$) implies the utility minimizing alter ego sets +Activating a concern about robustness, $\sigma < 0$, leads the utility-minimizing alter ego to set $$ -\tilde w_t = K(\sigma,\hat\beta)\,\mu_{s,t-1} -$$ +\tilde v_t = K(\sigma,\hat\beta)\,\mu_{s,t-1} +$$ (eq:rcs-K) -making the worst-case model for $\mu_{st}$: +making the worst-case model for $\mu_{st}$ $$ -\mu_{st} = (1 + \alpha\,K(\sigma,\hat\beta))\,\mu_{s,t-1} + \alpha\,\tilde\epsilon_t -$$ (eq:reversee3) +\mu_{st} = \zeta\,\mu_{s,t-1} + \alpha\,\tilde w_t, +\qquad \zeta \equiv 1 + \alpha\,K(\sigma,\hat\beta) +$$ (eq:rcs-scalar-worst) -For the allocation to remain the same, we require the robust Euler equation -$\hat\beta R\,\hat{\mathbb{E}}_t\mu_{s,t+1} = \mu_{st}$ to hold under the worst-case model, which gives +For the allocation to remain the same, we require the robust Euler equation $\hat\beta R\,\hat{\mathbb{E}}_t\mu_{s,t+1} = \mu_{st}$ to hold under the worst-case model, which gives $$ -(\hat\beta R)^{-1} = 1 + \alpha\, K(\sigma,\hat\beta) -$$ (eq:eulerdist) +\zeta = (\hat\beta R)^{-1} +$$ (eq:rcs-eulerdist) The minimizing agent's Bellman equation, a pure forecasting problem, yields $$ -\hat\zeta(\hat\beta) \equiv 1 + \alpha K(\sigma,\hat\beta) = \frac{1}{1 - \sigma\alpha^2 P(\hat\beta)} -$$ (eq:distort2) +\zeta = \frac{1}{1 - \sigma\alpha^2 P(\hat\beta)} +$$ (eq:rcs-zetaP) -where $P(\hat\beta)$ solves the scalar Bellman equation: +where $P(\hat\beta)$ solves the scalar Bellman equation $$ --P(\hat\beta) = \frac{\hat\beta - 1 + \sigma\alpha^2 + \sqrt{(\hat\beta-1+\sigma\alpha^2)^2 + 4\sigma\alpha^2}}{-2\sigma\alpha^2} -$$ (eq:distortcons) +P(\hat\beta) = \frac{\hat\beta - 1 + \sigma\alpha^2 + \sqrt{(\hat\beta-1+\sigma\alpha^2)^2 + 4\sigma\alpha^2}}{-2\sigma\alpha^2} +$$ (eq:rcs-riccati) -Solving {eq}`eq:eulerdist`-{eq}`eq:distortcons` for $\hat\beta$ gives exactly {eq}`eq:obseq`. +Solving {eq}`eq:rcs-eulerdist`-{eq}`eq:rcs-riccati` for $\hat\beta$ gives exactly {eq}`eq:rcs-oe`. ```` -Equation {eq}`eq:obseq` is the useful numerical object because it gives a straight-line map from the robustness parameter to the observationally equivalent discount factor. +Equation {eq}`eq:rcs-oe` is the useful numerical object because it gives a straight-line map from the robustness parameter to the observationally equivalent discount factor. ### Precautionary savings interpretation @@ -405,34 +428,27 @@ In the special case $\lambda = \delta_h = 0$, $s_t = c_t$ and the consumption ru $$ c_t = (1 - R^{-2}\beta^{-1})\!\left[Rk_{t-1} + \mathbb{E}_t\sum_{j=0}^{\infty}R^{-j}d_{t+j}\right] + \left(\frac{(R\beta)^{-1}-1}{R-1}\right)\!b -$$ (eq:consfunction) +$$ (eq:rcs-consfunction) -The **marginal propensity to consume** out of non-human wealth $Rk_{t-1}$ *equals* that out of -human wealth $\mathbb{E}_t\sum R^{-j}d_{t+j}$. +The **marginal propensity to consume** out of non-human wealth $Rk_{t-1}$ *equals* that out of human wealth $\mathbb{E}_t\sum R^{-j}d_{t+j}$. This equal-propensity property is a hallmark of the LQ model and *persists* when a concern for robustness is present, in contrast to usual precautionary-savings models with convex marginal utility. -{prf:ref}`thm-lqcs-oe1` says that with $\sigma < 0$, the observationally equivalent -$\hat\beta$ satisfies $\hat\beta < \beta$. +{prf:ref}`thm-rcs-oe1` says that with $\sigma < 0$, the observationally equivalent $\hat\beta$ satisfies $\hat\beta < \beta$. -If the starting point has $\beta R = 1$, then -$\hat\beta R < 1$. +If the starting point has $\beta R = 1$, then $\hat\beta R < 1$. -For a non-robust consumer with discount factor $\hat\beta$ at the same -interest rate, the Euler equation implies $\mathbb{E}_t c_{t+1} < c_t$: expected consumption -declines over time. +For a non-robust consumer with discount factor $\hat\beta$ at the same interest rate, the Euler equation implies $\mathbb{E}_t c_{t+1} < c_t$, so expected consumption declines over time. -This downward drift is the impatience offset in {prf:ref}`thm-lqcs-oe1`. +This downward drift is the impatience offset in {prf:ref}`thm-rcs-oe1`. It cancels the robust consumer's precautionary-savings motive, leaving the consumption and investment quantities unchanged. -The upward-drift comparison appears in {prf:ref}`thm-lqcs-oe2`, which asks the reverse observational-equivalence question. - -The classical precautionary motive arises because: +The classical precautionary motive arises because $$ u'''(c) > 0 \;\Rightarrow\; \mathbb{E}_t u'(c_{t+1}) > u'(\mathbb{E}_t c_{t+1}) \;\Rightarrow\; \mathbb{E}_t c_{t+1} > c_t -$$ +$$ (eq:rcs-prudence) This channel requires *convexity of marginal utility* and is absent with quadratic preferences. @@ -445,149 +461,55 @@ In contrast, the robustness-based precautionary motive operates through distorti The observational-equivalence result can be interpreted using a **Stackelberg multiplier game**. -After the minimizing agent has committed to a distortion process $\{w_{t+1}\}$, the maximizing consumer faces the following worst-case law of motion for the state $X_t$: +After the minimizing agent has committed to a distortion process $\{v_{t+1}\}$, the maximizing consumer faces the following worst-case law of motion for the state $X_t$: $$ \begin{aligned} -X_{t+1} &= \bigl(A - BF(\sigma,\hat\beta) + CK(\sigma,\hat\beta)\bigr) X_t + C\tilde\epsilon_{t+1} \\ +X_{t+1} &= \bigl(A - BF(\sigma,\hat\beta) + CK(\sigma,\hat\beta)\bigr) X_t + C\,w_{t+1} \\ \begin{pmatrix}b_t\\d_t\end{pmatrix} &= S X_t \end{aligned} -$$ (eq:sys2) +$$ (eq:rcs-worstcase-law) -A robust consumer with concerns about possible misspecification of the approximating model's stochastic process for non-financial income forms expectations of future income using the **distorted transition matrix** -$A - BF + CK$ rather than the approximating transition matrix $A - BF$. +A robust consumer forms expectations of future income using the **distorted transition matrix** $A - BF + CK$ rather than the approximating transition matrix $A - BF$. The distorted expectations operator $\hat{\mathbb{E}}_t$ satisfies $$ \hat{\mathbb{E}}_t X_{t+j} = (A - BF(\sigma,\hat\beta) + CK(\sigma,\hat\beta))^j X_t -$$ +$$ (eq:rcs-Ehat) Observational equivalence requires that the modified human-wealth formula $$ \hat\Psi_4 \sum_{j=0}^{\infty} R^{-j}\hat{\mathbb{E}}_t d_{t+j} -$$ +$$ (eq:rcs-humanwealth) equals its benchmark counterpart $\Psi_4 \sum_{j=0}^{\infty} R^{-j} \mathbb{E}_t d_{t+j}$. -This is achieved by a mutual adjustment of the coefficients $\hat\Psi_j$ through $\hat\beta$ and the distorted expectation operator $\hat{\mathbb{E}}_t$ through $\sigma$. +This is achieved by a mutual adjustment of the coefficients $\hat\Psi_j$ through $\hat\beta$ and of the distorted expectation operator $\hat{\mathbb{E}}_t$ through $\sigma$. The worst-case eigenvalue of $A - BF + CK$ exceeds that of $A - BF$ in modulus, so the worst-case distortions make the income process *more persistent* than under the approximating model. This is the precautionary motive in state-space form: the minimizing agent makes future income look more risky by introducing low-frequency persistence. -### Frequency domain interpretation - -```{index} single: Frequency Domain; permanent income model -``` - -The LQ permanent income framework has a natural frequency-domain interpretation. - -The consumer's concave utility makes him dislike **high-frequency** fluctuations in consumption, which he smooths by adjusting savings. - -High-frequency fluctuations are easier to smooth, so the consumer is automatically robust to misspecification of high-frequency features of the income process. - -**Low-frequency** fluctuations are harder to smooth because they are more persistent. - -In the frequency-domain notation of HST, the transfer function from shocks $\epsilon_t$ to the -target $s_t - b_t$ is $G(\zeta)$, and the frequency decomposition of the $H_2$ criterion is - -$$ -H_2 = -\frac{1}{2\pi}\int_{-\pi}^{\pi} \operatorname{trace}\!\bigl[G(\sqrt\beta\, e^{i\omega})^\top\,G(\sqrt\beta\, e^{i\omega})\bigr]\, d\omega -$$ - -The integrand $G^\top G$ is *largest at low frequencies* $\omega \approx 0$, where the consumer's welfare is most sensitive to income variability. - -Recognizing this, the minimizing agent concentrates the worst-case distortions at low frequencies. - -The distortion process has spectral density $W(\zeta)^\top W(\zeta)$ that is concentrated near $\omega = 0$. - -The variance of the worst-case shocks grows as $|\sigma|$ increases. - -### Detection error probabilities - -```{index} single: Detection Error Probabilities -``` - -A natural way to discipline the choice of $\sigma$ (or $\theta$) is to ask: **how difficult would -it be to statistically distinguish the approximating model from the worst-case model?** - -For a sample of length $T$, one can use a **log-likelihood ratio test** to compare the two -hypotheses. - -The **detection error probability** (DEP) is the probability of making the wrong -decision using the log-likelihood ratio statistic when one does not know which model generated the -data. - -Specifically: - -$$ -\text{DEP}(\sigma) = \frac{1}{2}\bigl[\mathbb{P}\{\text{prefer approx.} \mid \text{worst-case is true}\} - + \mathbb{P}\{\text{prefer worst-case} \mid \text{approx. is true}\}\bigr] -$$ - -When $\sigma = 0$ the two models are identical and DEP $= 0.5$. - -As $|\sigma|$ increases the -models diverge and the DEP falls toward zero. - -The full DEP calculation requires a specified approximating model, its worst-case counterpart, and the sample length used in the likelihood-ratio experiment. - -We compute such a DEP for a robust Bewley model below. - -```{note} -HST suggested that a DEP above 0.2 is "plausible", meaning the models are still hard enough to distinguish statistically that a concern for robustness is warranted. - -Values of $\sigma$ corresponding to DEP $\geq 0.2$ define a set of plausible worst-case models. -``` - -### Robustness of decision rules - -```{index} single: Robustness; payoff evaluation -``` - -To evaluate whether robust decision rules perform better than the non-robust rule when the data are -generated by a distorted model, define the **payoff** when the decision rule is designed for -robustness parameter $\sigma_2$ and the data are generated by the distorted model associated with -$\sigma_1$: - -$$ -\pi(\sigma_1;\sigma_2) = -\mathbb{E}_{0,\sigma_1}\sum_{t=0}^{\infty}\beta^t\, x_t^\top H(\sigma_2)^\top H(\sigma_2)\, x_t -$$ (eq:soln3) - -where the state evolves under decision rule $F(\sigma_2)$ and worst-case shocks $K(\sigma_1)$: - -$$ -x_{t+1} = \bigl(A - BF(\sigma_2) + CK(\sigma_1)\bigr)x_t + C\epsilon_{t+1} -$$ (eq:soln2) - -For $\sigma_1 = 0$ (approximating model generates data), the non-robust rule ($\sigma_2 = 0$) is -optimal by construction. - -As $\sigma_1$ decreases (the data are generated by increasingly -distorted models), the payoff of the $\sigma_2 = 0$ rule deteriorates faster than that of robust -rules. - -Computing the payoff comparison requires solving the full HST matrix problem for $F(\sigma_2)$ and $K(\sigma_1)$. +In the scalar reduction of the next section, this eigenvalue is $\zeta$, and we verify that $\zeta > 1$ while the approximating model has a unit root. ### Another observational equivalence result ```{index} single: Observational Equivalence; Theorem 2 ``` -````{prf:theorem} Observational Equivalence, II -:label: thm-lqcs-oe2 +````{prf:theorem} Observational equivalence, II +:label: thm-rcs-oe2 Fix all parameters except $(\sigma,\beta)$ and consider a consumption-investment allocation for $(\hat\sigma, \hat\beta)$ where $\hat\beta R = 1$ and $\hat\sigma < 0$. Then there exists $\tilde\beta > \hat\beta$ such that the $(\hat\sigma, \hat\beta)$ allocation also solves the $(0, \tilde\beta)$ problem. ```` -{prf:ref}`thm-lqcs-oe1` showed that starting from a benchmark with $\beta R = 1$, activating -robustness ($\sigma < 0$) is equivalent to *reducing* $\beta$. +{prf:ref}`thm-rcs-oe1` showed that starting from a benchmark with $\beta R = 1$, activating robustness is equivalent to *reducing* $\beta$. -{prf:ref}`thm-lqcs-oe2` goes in the opposite direction: it shows that the effects of activating a concern for robustness from a starting point with $\beta R = 1$ are replicated by *increasing* $\beta$ while setting $\sigma = 0$. +{prf:ref}`thm-rcs-oe2` goes in the opposite direction: the effects of activating a concern for robustness from a starting point with $\beta R = 1$ are replicated by *increasing* $\beta$ while setting $\sigma = 0$. In other words, when $\beta R = 1$, a concern for robustness operates like an *increase* in the discount factor, pushing $\beta R > 1$ and imparting an *upward drift* to the expected consumption profile. @@ -596,81 +518,69 @@ With $\hat\beta R = 1$ and $\hat\sigma < 0$, the robust Euler equation implies $$ \hat{\mathbb{E}}_t \mu_{c,t+1} = \mu_{ct} -$$ +$$ (eq:rcs-euler2) -One seeks $\tilde\beta > \hat\beta$ and $\sigma = 0$ such that the same allocation solves the -non-robust problem with discount factor $\tilde\beta$. +One seeks $\tilde\beta > \hat\beta$ and $\sigma = 0$ such that the same allocation solves the non-robust problem with discount factor $\tilde\beta$. -The key step is to observe that the worst-case distortion $K(\hat\sigma, \hat\beta)$ introduces a -drift in the marginal utility process that is equivalent to the drift produced by raising the -discount factor above $\hat\beta$. +The key step is to observe that the worst-case distortion $K(\hat\sigma, \hat\beta)$ introduces a drift in the marginal utility process that is equivalent to the drift produced by raising the discount factor above $\hat\beta$. Equating the two drifts and solving the scalar Bellman equation for $K$ yields $$ \tilde\beta(\hat\sigma) = \frac{\hat\beta(1+\hat\beta)}{2(1+\hat\sigma\alpha^2)} \left[1 + \sqrt{1 - 4\hat\beta\,\frac{1+\hat\sigma\alpha^2}{(1+\hat\beta)^2}}\right] -$$ (eq:obsequivn2) +$$ (eq:rcs-oe2) -The solution satisfies $\tilde\beta > \hat\beta$ when $\hat\sigma < 0$. +Setting $\hat\sigma = 0$ makes the square root equal $(1-\hat\beta)/(1+\hat\beta)$, so that $\tilde\beta = \hat\beta$. + +Since $1 + \hat\sigma\alpha^2$ decreases as $\hat\sigma$ falls below zero, both the prefactor and the square root increase, so $\tilde\beta > \hat\beta$ whenever $\hat\sigma < 0$. ```` -The map {eq}`eq:obsequivn2` is a closed form, so we can plot it directly. +### Comparing the two loci -The next figure compares the two observational-equivalence loci for the two-factor calibration, using $\alpha^2 = \sigma_1^2 + (1-\beta)^2\sigma_2^2$ (derived below in {eq}`eq:bew_alpha2`). +Both {eq}`eq:rcs-oe` and {eq}`eq:rcs-oe2` are closed forms, so we can plot them directly. -We start from a benchmark with $\hat\beta R = 1$, so $\hat\beta = \beta$. +We start from a benchmark with $\beta R = 1$. ```{code-cell} ipython3 --- mystnb: figure: caption: | - Two observational-equivalence experiments. Locus I (below $\beta$) holds the - *non-robust* agent fixed at $\beta R=1$ and reports the *robust* twin's - discount factor $\hat\beta(\sigma)$; locus II (above $\beta$) holds the - *robust* agent fixed at $\beta R=1$ and reports the *non-robust* twin's - discount factor $\tilde\beta(\hat\sigma)$. - name: fig-lqcs-oe-loci + Two observational-equivalence experiments. Locus I holds the *non-robust* + agent fixed at $\beta R = 1$ and reports the *robust* twin's discount + factor $\hat\beta(\sigma)$; locus II holds the *robust* agent fixed at + $\beta R = 1$ and reports the *non-robust* twin's discount factor + $\tilde\beta(\hat\sigma)$. + name: fig-rcs-oe-loci --- -β_bench = β # benchmark with β R = 1 -α2 = σ1**2 + (1 - β)**2 * σ2**2 # two-factor α^2 (see eq:bew_alpha2) - -σ_hat_vals = np.linspace(0.0, -0.16, 60) +σ_grid = np.linspace(0.0, -0.16, 60) -# Locus I (eq:obseq / eq:bew_locus): non-robust agent fixed at βR=1; -# report the robust twin's discount factor β̂(σ) < β -β_hat = β_bench + σ_hat_vals * α2 * β_bench / (1 - β_bench) +# locus I: non-robust agent fixed at βR=1, report the robust twin's β̂(σ) +β_hat = β + σ_grid * α2 * β / (1 - β) -# Locus II (eq:obsequivn2): robust agent fixed at βR=1; -# report the non-robust twin's discount factor β̃(σ̂) > β -disc = 1 - 4 * β_bench * (1 + σ_hat_vals * α2) / (1 + β_bench)**2 -β_tilde = (β_bench * (1 + β_bench)) / (2 * (1 + σ_hat_vals * α2)) \ - * (1 + np.sqrt(disc)) +# locus II: robust agent fixed at βR=1, report the non-robust twin's β̃(σ̂) +q = 1 + σ_grid * α2 +β_tilde = β * (1 + β) / (2 * q) * (1 + np.sqrt(1 - 4 * β * q / (1 + β)**2)) fig, ax = plt.subplots() -ax.plot(-σ_hat_vals, β_hat, lw=2, color='C3', - label=r'locus I: robust twin $\hat\beta(\sigma)<\beta$' - '\n(non-robust agent fixed at $\\beta R=1$)') -ax.plot(-σ_hat_vals, β_tilde, lw=2, color='C0', - label=r'locus II: non-robust twin $\tilde\beta(\hat\sigma)>\beta$' - '\n(robust agent fixed at $\\beta R=1$)') -ax.axhline(β_bench, color='k', linestyle=':', lw=1, +ax.plot(-σ_grid, β_hat, lw=2, color='C3', + label=r'locus I: robust twin $\hat\beta(\sigma) < \beta$') +ax.plot(-σ_grid, β_tilde, lw=2, color='C0', + label=r'locus II: non-robust twin $\tilde\beta(\hat\sigma) > \beta$') +ax.axhline(β, color='k', linestyle=':', lw=1, label=r'benchmark $\beta$ ($\beta R = 1$)') ax.set_xlabel(r'robustness concern $|\sigma|$') ax.set_ylabel('discount factor of the equivalent agent') -ax.legend(fontsize=8.5) +ax.legend() plt.show() - -print(f"at σ̂ = {σ_hat_vals[-1]:.3f}: β̃ = {β_tilde[-1]:.4f} > β = {β_bench}") -print(f" β̂ = {β_hat[-1]:.4f} < β = {β_bench}") ``` Both loci pass through the benchmark $\beta$ at $\sigma = 0$ and separate as the robustness concern grows. -The key to reading the figure is that the two loci hold *different* agents fixed, so the discount factor plotted on the vertical axis refers to a different agent on each curve. +The key to reading {numref}`fig-rcs-oe-loci` is that the two loci hold *different* agents fixed, so the discount factor plotted on the vertical axis refers to a different agent on each curve. -Locus I, from {prf:ref}`thm-lqcs-oe1`, holds the **non-robust** agent fixed at the benchmark $(\sigma = 0, \beta)$ with $\beta R = 1$ and reports the discount factor $\hat\beta(\sigma) < \beta$ of the **robust** agent that mimics it. +Locus I, from {prf:ref}`thm-rcs-oe1`, holds the **non-robust** agent fixed at the benchmark $(\sigma = 0, \beta)$ with $\beta R = 1$ and reports the discount factor $\hat\beta(\sigma) < \beta$ of the **robust** agent that mimics it. This is the sense in which HST call a concern for robustness observationally equivalent to a *lower* discount factor: because robustness already makes the agent save more, its discount factor must be lowered to hold the allocation at the benchmark. @@ -680,488 +590,640 @@ The robust twin chooses the identical consumption process, so it too satisfies $ The lower $\hat\beta$, which has $\hat\beta R < 1$, would on its own impart a downward drift, but the robust agent's precautionary saving offsets it exactly, leaving expected consumption flat. -Locus II, from {prf:ref}`thm-lqcs-oe2`, instead holds the **robust** agent fixed at $(\hat\sigma, \beta)$ with $\beta R = 1$ and reports the discount factor $\tilde\beta(\hat\sigma) > \beta$ of the **non-robust** agent that mimics it. +Locus II, from {prf:ref}`thm-rcs-oe2`, instead holds the **robust** agent fixed at $(\hat\sigma, \beta)$ with $\beta R = 1$ and reports the discount factor $\tilde\beta(\hat\sigma) > \beta$ of the **non-robust** agent that mimics it. Here there is no impatience offset, so the common allocation inherits the robust agent's precautionary *upward* drift, which the non-robust twin reproduces through $\tilde\beta R > 1$. The two experiments encode the *same* economics: a concern for robustness adds precautionary saving that acts like extra patience. -They differ only in which agent is anchored at $\beta R = 1$, and hence in whether the common saving motive shows up as an exactly-offsetting impatience adjustment (locus I, expected consumption flat) or as an upward drift in expected consumption (locus II). +They differ only in which agent is anchored at $\beta R = 1$, and hence in whether the common saving motive shows up as an exactly-offsetting impatience adjustment, as in locus I where expected consumption is flat, or as an upward drift in expected consumption, as in locus II. -### A robust LQ Bewley model +(rcs-scalar)= +## The scalar worst-case model -```{index} single: Robust Bewley Model +```{index} single: Robust Control; breakdown point ``` -We now synthesise the lecture by embedding the Bewley economy of {doc}`lq_bewley_complete_markets` into the HST framework and applying the observational-equivalence theorem. +The proof of {prf:ref}`thm-rcs-oe1` reduced the robust problem to a scalar forecasting problem in the marginal utility $\mu_{st}$. -We shall construct a family of **robust Bewley economies**, parameterised by a robustness level $\sigma \leq 0$, whose equilibrium quantities are identical to those of the plain vanilla Bewley model. +That scalar problem can be solved in closed form, which is worth doing because it makes the worst-case dynamics, the breakdown point, and the frequency-domain results all transparent. -We first map the Bewley economy into HST notation, specialising the robust model to -$\lambda = \delta_h = 0$ (no habits, no durable goods) and to a -pure endowment economy (no physical capital, $k_t = 0$). +### A closed-form solution -In this case: +Combining {eq}`eq:rcs-eulerdist` with {eq}`eq:rcs-oe` and $R = \beta^{-1}$ gives the worst-case persistence directly: -Services equal consumption: $s_t = c_t$. +$$ +\zeta(\sigma) = \frac{1}{\hat\beta(\sigma) R} = \frac{\beta}{\hat\beta(\sigma)} += \left[1 + \frac{\sigma\alpha^2}{1-\beta}\right]^{-1} +$$ (eq:rcs-zeta) -The only traded security is the one-period risk-free bond, and we write the household's net asset position as $a_t=-b_t$ so that positive $a_t$ denotes wealth rather than debt. +This is the central formula of the lecture. -The endowment process follows the state-space representation {eq}`eq:rs-endowment`. +**The worst-case persistence of marginal utility is exactly the ratio of the two discount factors.** -The household's augmented state vector is $x_t = [a_t,\; z_t^\top]^\top$, and the law of motion -{eq}`eq:law0` specialises to +Since $\hat\beta(\sigma) < \beta$ for $\sigma<0$, we have $\zeta(\sigma) > 1$: the approximating model for $\mu_{st}$ has a unit root, and the worst-case model is mildly explosive. -$$ -\begin{pmatrix} a_{t+1} \\ z_{t+1} \end{pmatrix} -= -\underbrace{\begin{pmatrix} R & R\check{G} \\ 0 & \check{A} \end{pmatrix}}_{A} -\begin{pmatrix} a_t \\ z_t \end{pmatrix} -+ -\underbrace{\begin{pmatrix} -R \\ 0 \end{pmatrix}}_{B} -c_t -+ -\underbrace{\begin{pmatrix} 0 \\ \check{C} \end{pmatrix}}_{C} -\epsilon_{t+1} -$$ (eq:bew_law) +That is the scalar counterpart of the statement in {eq}`eq:rcs-worstcase-law` that the worst-case transition matrix has a larger eigenvalue than the approximating one. -The objective is $\mathbb{E}_0 \sum_{t=0}^\infty \beta^t [-(c_t - \gamma)^2/2]$, which is the HST -criterion {eq}`eq:income5` with $\sigma = 0$ and $b_t \equiv \gamma$ (a fixed bliss level). +We can also solve {eq}`eq:rcs-riccati` explicitly. -The robust Bellman equation {eq}`eq:income1` with $\sigma = 0$ therefore reduces exactly to -the LQ problem of {doc}`lq_permanent_income`, confirming that the HST framework nests the Bewley model. +Write $u = \sigma\alpha^2$ and $\delta = 1-\beta$, so that {eq}`eq:rcs-oe` reads $\hat\beta - 1 + u = (u - \delta^2)/\delta$. -We next compute the robustness parameter $\alpha^2$. +The discriminant in {eq}`eq:rcs-riccati` is then a **perfect square**: -From the $(c_t,z_t)$ representation {eq}`eq:rs-crep`, the consumption innovation is +$$ +(\hat\beta - 1 + u)^2 + 4u = \frac{(u-\delta^2)^2}{\delta^2} + 4u = \left(\frac{u+\delta^2}{\delta}\right)^2 +$$ (eq:rcs-disc) + +so the two roots of {eq}`eq:rcs-riccati` are available in closed form: $$ -c_{t+1} - c_t = h\, w_{t+1}, \qquad -h = (1-\beta)\,\check{G}(I-\beta\check{A})^{-1}\check{C} -$$ (eq:bew_cinno) +P = -\frac{1}{1-\beta} +\qquad\text{and}\qquad +P = \frac{1-\beta}{\sigma\alpha^2} +$$ (eq:rcs-roots) -The vector $h$ plays the role of $\nu^\top = M_s C$ in the HST scalar $\alpha$ formula -{eq}`eq:hsoffset2`. +Substituting into {eq}`eq:rcs-zetaP`, the first root reproduces {eq}`eq:rcs-zeta` while the second gives the constant $\zeta = R$, which violates the Euler equation {eq}`eq:rcs-eulerdist` except at the single point where the two roots coincide. -Consequently, +So the economically relevant root is the constant $P = -(1-\beta)^{-1}$, independent of $\sigma$. -$$ -\alpha^2 = h h^\top = (1-\beta)^2\, -\check{G}(I-\beta\check{A})^{-1}\check{C}\check{C}^\top(I-\beta\check{A}^\top)^{-1}\check{G}^\top -$$ (eq:bew_alpha) +```{note} +Selecting the root numerically, by taking whichever of the two is closer to the target $(\hat\beta R)^{-1}$, is treacherous. -For the two-factor model {eq}`eq:pi-twofactor` with $\check{A} = \mathrm{diag}(1,0)$ and -$\check{C} = \mathrm{diag}(\sigma_1,\sigma_2)$ this simplifies to +Because {eq}`eq:rcs-disc` is a perfect square, the square root is $|u+\delta^2|/\delta$, and the *sign* in front of it that selects $P = -(1-\beta)^{-1}$ flips as $u$ crosses $-\delta^2$. -$$ -\alpha^2 = \sigma_1^2 + (1-\beta)^2\,\sigma_2^2 -$$ (eq:bew_alpha2) +A solver that picks the root by a distance criterion will silently switch branches there. +``` + +### The breakdown point + +{prf:ref}`thm-rcs-oe1` asserts the existence of a lower bound $\underline\sigma < 0$ without saying what it is. -The permanent shock variance $\sigma_1^2$ enters with coefficient 1 because a unit permanent -shock is *fully* capitalised into consumption. +The scalar model tells us. -The transitory shock variance $\sigma_2^2$ -enters with the small coefficient $(1-\beta)^2$ because only its annuity value is consumed. +The minimizing agent's problem has a finite value only if the discounted worst-case state is square summable, that is, only if $\hat\beta\,\zeta^2 < 1$. -Applying {prf:ref}`thm-lqcs-oe1` {eq}`eq:obseq` with equilibrium interest rate $R = \beta_0^{-1}$ and -$\alpha^2$ from {eq}`eq:bew_alpha2` gives the **Bewley observational equivalence locus**: +Using $\hat\beta = \beta/\zeta$ from {eq}`eq:rcs-zeta`, this condition is $\beta\zeta < 1$, or equivalently $\zeta < R$. + +Substituting {eq}`eq:rcs-zeta` and solving gives the **breakdown point** $$ -\hat\beta(\sigma) = \beta_0 + \frac{\sigma\,\alpha^2\,\beta_0}{1-\beta_0} -$$ (eq:bew_locus) +\underline\sigma = -\frac{(1-\beta)^2}{\alpha^2} +$$ (eq:rcs-breakdown) -For $\sigma < 0$, we have $\hat\beta(\sigma) < \beta_0$. +At $\sigma = \underline\sigma$ three things happen at once. -An agent with the pair -$(\sigma, \hat\beta(\sigma))$ on this locus is more concerned about model misspecification -(lower $\sigma$) but also more impatient (lower $\hat\beta$); the two forces cancel exactly, -leaving the consumption decision rule unchanged. +- The discriminant {eq}`eq:rcs-disc` has a double root, so the two roots in {eq}`eq:rcs-roots` coincide. +- The worst-case persistence reaches $\zeta = R$, so $\hat\beta\zeta^2 = 1$ exactly. +- The observationally equivalent discount factor reaches $\hat\beta = \beta^2$. -These ingredients combine into a robust Bewley equilibrium. +For $\sigma < \underline\sigma$ the robust control problem has no solution, and any numerical answer reported there is meaningless. -````{prf:proposition} -:label: prop-lqcs-bewley +We therefore restrict every plot below to $\sigma \in (\underline\sigma, 0]$. -Suppose all agents in the Bewley economy share a common pair -$(\sigma, \hat\beta(\sigma))$ lying on the locus {eq}`eq:bew_locus`, with $R = \beta_0^{-1}$. +### Verifying the closed form -Then every agent's optimal consumption plan is identical to that of the plain vanilla -$(\sigma = 0,\, \beta_0)$ economy, and the equilibrium interest rate remains $R = \beta_0^{-1}$. -```` +The following cell solves the quadratic {eq}`eq:rcs-riccati` numerically and checks it against the closed forms {eq}`eq:rcs-zeta` and {eq}`eq:rcs-roots`. -````{prf:proof} -By {prf:ref}`thm-lqcs-oe1`, each agent's consumption-saving rule is identical to the benchmark. +```{code-cell} ipython3 +def worst_case_persistence(σ, β, α2): + """ + Worst-case persistence ζ(σ) of marginal utility on the + observational-equivalence locus, from eq:rcs-zeta. + """ + return 1 / (1 + σ * α2 / (1 - β)) -The goods-market clearing condition $\int c_t^i\, di = Y$ is therefore satisfied at -$R = \beta_0^{-1}$ for the same reason as in the benchmark Bewley economy. -```` -#### Heterogeneous $(\beta_i, \sigma_i)$ preferences +def solve_scalar_riccati(σ, β, α2): + """ + Solve the scalar Bellman equation eq:rcs-riccati by brute force and + return both roots together with the implied persistence ζ = 1/(1-σα²P). + """ + β_hat = β + σ * α2 * β / (1 - β) + u = σ * α2 + disc = (β_hat - 1 + u)**2 + 4 * u + roots = np.array([(β_hat - 1 + u + s * np.sqrt(disc)) / (-2 * u) + for s in (1.0, -1.0)]) + return roots, 1 / (1 - u * roots) + + +σ_lo = -(1 - β)**2 / α2 # breakdown point, eq:rcs-breakdown +print(f"breakdown point σ̲ = {σ_lo:.6f}") +print(f"there β̂ = {β + σ_lo * α2 * β / (1 - β):.6f} " + f"(β² = {β**2:.6f})") +print(f" ζ = {worst_case_persistence(σ_lo, β, α2):.6f} " + f"(R = {R:.6f})") + +print(f"\n{'σ':>10}{'P (numerical roots)':>28}{'-1/(1-β)':>12}" + f"{'ζ (num)':>12}{'ζ (closed)':>12}") +for σ in [-0.02, -0.05, -0.09, -0.105]: + roots, ζs = solve_scalar_riccati(σ, β, α2) + keep = np.argmin(np.abs(roots + 1 / (1 - β))) + print(f"{σ:10.3f}{str(np.round(roots, 4)):>28}{-1 / (1 - β):12.4f}" + f"{ζs[keep]:12.6f}{worst_case_persistence(σ, β, α2):12.6f}") +``` + +One root is pinned at $-(1-\beta)^{-1} = -20$ for every $\sigma$, exactly as {eq}`eq:rcs-roots` predicts, and the implied $\zeta$ agrees with the closed form to displayed precision. + +Notice also that the two roots approach each other as $\sigma$ falls toward $\underline\sigma \approx -0.11$. -A richer extension populates the economy with a **continuum of types**, each indexed by a -robustness parameter $\sigma_i \in [\underline\sigma, 0]$, with discount factor +The next figure plots the worst-case impulse response $\zeta^h$ over the admissible range. + +```{code-cell} ipython3 +--- +mystnb: + figure: + caption: | + Worst-case impulse response of marginal utility. The approximating model + has a unit root; the worst-case model is increasingly explosive as + $\sigma$ falls toward the breakdown point $\underline\sigma$. + name: fig-rcs-irf +--- +horizons = np.arange(31) + +fig, ax = plt.subplots() +for frac in [0.0, 0.4, 0.8, 0.98]: + σ = frac * σ_lo + ζ = worst_case_persistence(σ, β, α2) + ax.plot(horizons, ζ**horizons, lw=2, + label=rf'$\sigma={σ:.3f}$, $\zeta={ζ:.3f}$') + +ax.set_xlabel('horizon $h$') +ax.set_ylabel(r'response of $\mu_{s,t+h}$') +ax.legend() +plt.show() +``` + +{numref}`fig-rcs-irf` shows that as $\sigma$ falls, the minimizing agent converts the unit root of the approximating model into mild explosiveness. + +At the breakdown point the response would grow at the gross interest rate $R$, so that the discounted worst-case state ceases to be square summable. + +## Frequency domain interpretation + +```{index} single: Frequency Domain; permanent income model +``` + +The LQ permanent income framework has a natural frequency-domain interpretation. + +The consumer's concave utility makes him dislike **high-frequency** fluctuations in consumption, which he smooths by adjusting savings. + +High-frequency fluctuations are easy to smooth, so the consumer is automatically robust to misspecification of high-frequency features of the income process. + +**Low-frequency** fluctuations are harder to smooth because they are more persistent. + +In the frequency-domain notation of HST, the transfer function from shocks $w_t$ to the target $s_t - b_t$ is $G(\cdot)$, and the frequency decomposition of the $H_2$ criterion is $$ -\beta_i = \hat\beta(\sigma_i) = \beta_0 + \frac{\sigma_i\,\alpha^2\,\beta_0}{1-\beta_0} -$$ (eq:bew_heterog) +H_2 = -\frac{1}{2\pi}\int_{-\pi}^{\pi} \operatorname{trace}\!\bigl[G(\sqrt\beta\, e^{i\omega})^\top\,G(\sqrt\beta\, e^{i\omega})\bigr]\, d\omega +$$ (eq:rcs-h2) -Since all pairs $(\sigma_i, \beta_i)$ lie on {eq}`eq:bew_locus`, every agent adopts the **same consumption rule** as the benchmark. +The evaluation at $\sqrt{\beta}\,e^{i\omega}$ rather than at $e^{i\omega}$ is essential and not merely conventional. -Aggregate dynamics are unchanged because the cross-section mean of consumption equals $Y$ and the cross-section variance grows at rate $\alpha^2$ per period. +Both the approximating model and the worst-case model for $\mu_{st}$ are non-stationary, the first with a unit root and the second explosive, so neither has an ordinary spectral density. -The equilibrium interest rate is unchanged: $R = \beta_0^{-1}$. +Discounting by $\sqrt\beta$ is exactly what makes the object in {eq}`eq:rcs-h2` finite. -Agents are observationally indistinguishable to an outside econometrician because data on $(c_t^i, a_t^i)$ cannot reveal whether agent $i$ has $\sigma_i = 0$ or $\sigma_i < 0$. +In the scalar model of {ref}`rcs-scalar`, the target is $s_t - b_t = -\mu_{st}$, and from {eq}`eq:rcs-scalar-worst` we can read off the transfer function and the discounted spectral density -Agents differ in their internal model because an agent with $\sigma_i < 0$ applies a worst-case distortion $w_{t+1}^i = K(\sigma_i, \beta_i)\,\mu_{s,t}^i$ to her conditional expectations, while an agent with $\sigma_i = 0$ takes the approximating model at face value. +$$ +G(z) = \frac{\alpha}{1-\zeta z}, +\qquad +S(\omega;\sigma) = \bigl|G(\sqrt{\hat\beta}\, e^{i\omega})\bigr|^2 += \frac{\alpha^2}{\bigl|1 - \zeta\sqrt{\hat\beta}\, e^{i\omega}\bigr|^2} +$$ (eq:rcs-spectrum) + +This is finite precisely when $\zeta\sqrt{\hat\beta} < 1$, which is the condition $\hat\beta\zeta^2<1$ that defines the breakdown point {eq}`eq:rcs-breakdown`. -This sets the stage for a Bewley model with **heterogeneous ambiguity aversion**: although -every agent acts identically in terms of observable choices, they hold different subjective -models of the income process and have different attitudes toward model uncertainty. +So the frequency-domain object and the breakdown point are two views of the same restriction. -#### Computation +At $\sigma = 0$ we have $\zeta = 1$ and $\hat\beta = \beta$, and {eq}`eq:rcs-spectrum` reduces to the discounted spectral density of a random walk. ```{code-cell} ipython3 -# Bewley parameters -β0_bew = β # 0.95 -σ1_bew = σ1 # 0.15 -σ2_bew = σ2 # 0.30 -R_bew = 1.0 / β0_bew - -# Two-factor Bewley α^2 -α2_bew = σ1_bew**2 + (1 - β0_bew)**2 * σ2_bew**2 - -print(f"α^2 (Bewley, two-factor) = {α2_bew:.6f}") -print(f" permanent component σ1^2 = {σ1_bew**2:.6f} " - f"({100*σ1_bew**2/α2_bew:.1f} % of α^2)") -print(f" transitory component (1-β)^2σ2^2= {(1-β0_bew)**2*σ2_bew**2:.6f} " - f"({100*(1-β0_bew)**2*σ2_bew**2/α2_bew:.1f} % of α^2)") +--- +mystnb: + figure: + caption: | + Left: the discounted spectral density of the robust consumer's target. + Right: the same densities relative to the approximating model. The + minimizing agent concentrates its distortion at low frequencies, where + the permanent income consumer is least able to smooth. + name: fig-rcs-spectrum +--- +ω = np.linspace(0, np.pi, 400) + + +def spectrum(σ, β, α2): + ζ = worst_case_persistence(σ, β, α2) + β_hat = β / ζ + return α2 / np.abs(1 - ζ * np.sqrt(β_hat) * np.exp(1j * ω))**2 + + +S0 = spectrum(0.0, β, α2) + +fig, axes = plt.subplots(1, 2, figsize=(11, 4)) + +for frac in [0.0, 0.4, 0.8, 0.95]: + σ = frac * σ_lo + S = spectrum(σ, β, α2) + axes[0].plot(ω, S, lw=2, label=rf'$\sigma={σ:.3f}$') + axes[1].plot(ω, S / S0, lw=2, label=rf'$\sigma={σ:.3f}$') + +axes[0].set_yscale('log') +axes[0].set_xlabel(r'frequency $\omega$') +axes[0].set_ylabel(r'$S(\omega;\sigma)$') +axes[0].set_title('discounted spectral density') +axes[0].legend() + +axes[1].set_yscale('log') +axes[1].set_xlabel(r'frequency $\omega$') +axes[1].set_ylabel(r'$S(\omega;\sigma)\,/\,S(\omega;0)$') +axes[1].set_title('relative to the approximating model') +axes[1].legend() + +fig.tight_layout() +plt.show() ``` -The calculation shows why permanent shocks dominate $\alpha^2$ in this calibration. +The left panel of {numref}`fig-rcs-spectrum` shows that $S(\omega;\sigma)$ is largest at $\omega \approx 0$, where the consumer's welfare is most sensitive to income variability. -We now solve the scalar robust forecasting problem attached to this $\alpha^2$. +The right panel makes the key point: the ratio $S(\omega;\sigma)/S(\omega;0)$ is sharply decreasing in $\omega$. -The solution selects the Bellman-equation root that satisfies the observational-equivalence Euler equation. +Recognizing where the consumer is vulnerable, the minimizing agent concentrates the worst-case distortions at low frequencies, and it does so more aggressively as $|\sigma|$ grows. -```{code-cell} ipython3 -def robust_scalar_solution(σ, β0, α2): - """ - Solve the scalar robust marginal-utility problem on the - observational-equivalence locus. - """ - α = np.sqrt(α2) - R = 1.0 / β0 +The peak at $\omega = 0$ equals $\alpha^2/(1-\sqrt{\beta\zeta})^2$ and diverges as $\sigma \to \underline\sigma$, which is another way of seeing the breakdown point. - if np.isclose(σ, 0.0): - return β0, np.nan, 1.0, 0.0 +Because the distortion is $v_t = K\mu_{s,t-1}$ by {eq}`eq:rcs-K`, the spectral density of the distortion process is $K^2$ times the density of $\mu_{s,t-1}$, so it inherits the same low-frequency concentration. - β_hat = β0 + σ * α2 * β0 / (1 - β0) - disc = (β_hat - 1 + σ * α2)**2 + 4 * σ * α2 - root_disc = np.sqrt(max(disc, 0.0)) - target_ζ = 1 / (β_hat * R) +## Detection error probabilities - candidates = [] - for sign in (1.0, -1.0): - P = (β_hat - 1 + σ * α2 + sign * root_disc) / (-2 * σ * α2) - ζ = 1 / (1 - σ * α2 * P) - K = (ζ - 1) / α - candidates.append((abs(ζ - target_ζ), P, ζ, K)) +```{index} single: Detection Error Probabilities +``` - _, P, ζ, K = min(candidates, key=lambda x: x[0]) - return β_hat, P, ζ, K +A natural way to discipline the choice of $\sigma$ is to ask: **how difficult would it be to distinguish the approximating model from the worst-case model statistically?** +For a sample of length $T$, one can use a **log-likelihood ratio test** to compare the two hypotheses. -def log_likelihood_ratio(paths, ζ, α): - """ - Return log p_worst(path) - log p_approx(path). - """ - lag = paths[:, :-1] - lead = paths[:, 1:] - ll_worst = -0.5 * np.sum(((lead - ζ * lag) / α)**2, axis=1) - ll_approx = -0.5 * np.sum(((lead - lag) / α)**2, axis=1) - return ll_worst - ll_approx +The **detection error probability** (DEP) is the probability of making the wrong decision using the log-likelihood ratio statistic when one does not know which model generated the data: + +$$ +\mathrm{DEP}(\sigma) = \frac{1}{2}\bigl[\mathbb{P}\{\text{prefer approximating} \mid \text{worst-case is true}\} + + \mathbb{P}\{\text{prefer worst-case} \mid \text{approximating is true}\}\bigr] +$$ (eq:rcs-dep) + +When $\sigma = 0$ the two models are identical and $\mathrm{DEP} = 0.5$. + +As $|\sigma|$ increases the models diverge and the DEP falls toward zero. +In the scalar model the two hypotheses are fully explicit: -def simulate_scalar_paths(ζ, α, T, n_paths, seed): +$$ +\text{approximating:}\quad \mu_{t+1} = \mu_t + \alpha w_{t+1}, +\qquad +\text{worst-case:}\quad \mu_{t+1} = \zeta(\sigma)\,\mu_t + \alpha w_{t+1} +$$ (eq:rcs-hypotheses) + +Both have Gaussian innovations with the same variance $\alpha^2$, so the log-likelihood ratio is a difference of sums of squares. + +```{code-cell} ipython3 +def simulate_paths(ζ, α, T, n_paths, seed): + "Simulate n_paths draws of μ_{t+1} = ζ μ_t + α w_{t+1} from μ_0 = 0." rng = np.random.default_rng(seed) paths = np.zeros((n_paths, T + 1)) shocks = rng.standard_normal((n_paths, T)) - for t in range(T): paths[:, t + 1] = ζ * paths[:, t] + α * shocks[:, t] - return paths +def log_likelihood_ratio(paths, ζ, α): + "Return log p_worst(path) - log p_approx(path)." + lag, lead = paths[:, :-1], paths[:, 1:] + return 0.5 * (np.sum(((lead - lag) / α)**2, axis=1) + - np.sum(((lead - ζ * lag) / α)**2, axis=1)) + + def detection_error_probability(ζ, α, T=40, n_paths=10_000, seed=1234): - """ - Finite-sample DEP for the approximating and worst-case scalar laws. - """ + "Finite-sample DEP for the approximating and worst-case scalar laws." if np.isclose(ζ, 1.0): return 0.5 + approx = simulate_paths(1.0, α, T, n_paths, seed) + worst = simulate_paths(ζ, α, T, n_paths, seed + 1) + return 0.5 * (np.mean(log_likelihood_ratio(worst, ζ, α) < 0) + + np.mean(log_likelihood_ratio(approx, ζ, α) > 0)) +``` - approx_paths = simulate_scalar_paths(1.0, α, T, n_paths, seed) - worst_paths = simulate_scalar_paths(ζ, α, T, n_paths, seed + 1) +We use $T = 40$, which is ten years of quarterly data. - llr_approx = log_likelihood_ratio(approx_paths, ζ, α) - llr_worst = log_likelihood_ratio(worst_paths, ζ, α) +Before plotting, it is worth noting a property that makes the DEP the right way to report a robustness concern. - return 0.5 * (np.mean(llr_worst < 0) + np.mean(llr_approx > 0)) +The parameter $\sigma$ is *not* scale free: it always enters through the product $\sigma\alpha^2$, and $\alpha^2$ has the units of consumption squared. + +Doubling the units in which consumption is measured therefore changes the numerical value of $\sigma$ that represents a given concern for robustness. + +The DEP has no such problem, because $\alpha$ cancels from the likelihood ratio in {eq}`eq:rcs-hypotheses`. + +```{code-cell} ipython3 +ζ_test = worst_case_persistence(0.6 * σ_lo, β, α2) +for scale in [0.5, 1.0, 2.0]: + dep = detection_error_probability(ζ_test, scale * α) + print(f"α scaled by {scale:>4}: DEP = {dep:.4f}") ``` -The next figure reports worst-case dynamics and model-detection probabilities implied by this solved scalar problem. +The DEP depends only on $\zeta$ and on the sample size $T$. ```{code-cell} ipython3 --- mystnb: figure: - caption: Solved robust scalar model - name: fig-lqcs-robust-scalar + caption: | + Finite-sample detection error probability against the robustness + concern, over the admissible range $(\underline\sigma, 0]$, for two + sample lengths. Values below the dashed line are the ones HST regard as + implausible, because the approximating and worst-case models would then + be too easy to tell apart. + name: fig-rcs-dep --- -α_bew = np.sqrt(α2_bew) -β_min = 0.88 -σ_min = (β_min - β0_bew) * (1 - β0_bew) / (α2_bew * β0_bew) -σ_vals = np.linspace(0.0, σ_min, 31) - -solutions = np.array([robust_scalar_solution(σ, β0_bew, α2_bew) for σ in σ_vals]) -β_hat_vals = solutions[:, 0] -ζ_vals = solutions[:, 2] -K_vals = solutions[:, 3] -dep_vals = np.array([ - detection_error_probability(ζ, α_bew) - for ζ in ζ_vals -]) - -fig, axes = plt.subplots(1, 2, figsize=(11.2, 4.1)) +σ_vals = np.linspace(0.0, 0.999 * σ_lo, 31) +ζ_vals = worst_case_persistence(σ_vals, β, α2) -horizons = np.arange(31) -for σ in [0.0, σ_vals[10], σ_vals[20]]: - β_hat, P, ζ, K = robust_scalar_solution(σ, β0_bew, α2_bew) - label = rf'$\sigma={σ:.3f}$, $\zeta={ζ:.3f}$' - axes[0].plot(horizons, ζ**horizons, lw=2, label=label) - -axes[0].set_xlabel('horizon') -axes[0].set_ylabel(r'response of $\mu_{s,t+h}$') -axes[0].set_title('worst-case impulse response') -axes[0].legend(fontsize=8.5) - -axes[1].plot(-σ_vals, dep_vals, lw=2, color='C0') -axes[1].axhline(0.2, color='C3', linestyle='--', lw=1.2, - label='DEP = 0.2') -axes[1].set_xlabel(r'robustness concern $-\sigma$') -axes[1].set_ylabel('detection error probability') -axes[1].set_ylim(0.0, 0.52) -axes[1].set_title('finite-sample detectability') -axes[1].legend(fontsize=8.5) - -fig.tight_layout() +fig, ax = plt.subplots() +for T, color in [(40, 'C0'), (160, 'C2')]: + dep_vals = np.array([detection_error_probability(ζ, α, T=T) + for ζ in ζ_vals]) + ax.plot(-σ_vals, dep_vals, lw=2, color=color, label=f'$T = {T}$') + +ax.axhline(0.2, color='C3', linestyle='--', lw=1.2, label='DEP = 0.2') +ax.axvline(-σ_lo, color='k', linestyle=':', lw=1, + label='breakdown point') +ax.set_xlabel(r'robustness concern $|\sigma|$') +ax.set_ylabel('detection error probability') +ax.set_ylim(0.0, 0.52) +ax.legend() plt.show() ``` -The left panel shows that the solved worst-case law makes marginal utility more persistent as $\sigma$ becomes more negative. +```{note} +HST suggested that a DEP above 0.2 is "plausible", meaning the models are still hard enough to distinguish statistically that a concern for robustness is warranted. + +Values of $\sigma$ with $\mathrm{DEP} \geq 0.2$ therefore define a set of plausible worst-case models. +``` -The right panel computes the DEP from the exact likelihood ratio between the approximating scalar law $\mu_{t+1}=\mu_t+\alpha\epsilon_{t+1}$ and the solved worst-case law $\mu_{t+1}=\zeta(\sigma)\mu_t+\alpha\epsilon_{t+1}$. +{numref}`fig-rcs-dep` shows that the two disciplines interact in an interesting way. -### Concluding remarks +At $T = 40$ the DEP curve reaches the $0.2$ threshold essentially exactly at the breakdown point. + +In other words, with ten years of quarterly data, every robustness concern that this model can represent at all is also statistically plausible: the mathematical limit binds first, and the statistical one is slack. + +```{code-cell} ipython3 +for frac in [0.6, 0.9, 0.999]: + ζ = worst_case_persistence(frac * σ_lo, β, α2) + print(f"σ = {frac * σ_lo:.5f} (= {frac:.3f} σ̲): " + f"DEP = {detection_error_probability(ζ, α):.4f}") +``` + +That near-coincidence is a property of this calibration and of $T = 40$, not a theorem. + +At $T = 160$ the ordering reverses, and statistical detectability binds well before the breakdown point. + +```{code-cell} ipython3 +σ_report = 0.6 * σ_lo +ζ_report = worst_case_persistence(σ_report, β, α2) +print(f"at σ = {σ_report:.4f} (ζ = {ζ_report:.4f}):") +for T in [20, 40, 80, 160]: + print(f" T = {T:>3}: DEP = " + f"{detection_error_probability(ζ_report, α, T=T):.4f}") +``` + +The lesson is that the plausible amount of model uncertainty depends on how much data the econometrician is imagined to have. + +## Robustness of decision rules + +```{index} single: Robustness; payoff evaluation +``` + +A natural follow-up question is whether robust decision rules perform better than the non-robust rule when the data are in fact generated by a distorted model. + +Define the **payoff** when the decision rule is designed for robustness parameter $\sigma^r$ and the data are generated by the distorted model associated with $\sigma^d$: + +$$ +\pi(\sigma^d;\sigma^r) = -\mathbb{E}_{0,\sigma^d}\sum_{t=0}^{\infty}\beta^t\, x_t^\top H(\sigma^r)^\top H(\sigma^r)\, x_t +$$ (eq:rcs-payoff) + +where the state evolves under decision rule $F(\sigma^r)$ and worst-case shocks $K(\sigma^d)$: + +$$ +x_{t+1} = \bigl(A - BF(\sigma^r) + CK(\sigma^d)\bigr)x_t + C\,w_{t+1} +$$ (eq:rcs-payoff-law) -We close with a summary of the key messages from all three lectures. +There is an important caveat about *which* family of rules to compare, and it follows directly from {prf:ref}`thm-rcs-oe1`. -The LQ permanent income model, a rational-expectations version of Friedman's permanent income hypothesis, has two complementary state-space representations: +Along the observational-equivalence locus $(\sigma, \hat\beta(\sigma))$, every agent chooses the *same* decision rule. -1. **$(b_t, z_t)$ representation**: emphasises that the consumer's optimal borrowing is history - dependent and cointegrated with consumption. +So $F(\sigma^r)$ does not vary along the locus, and $\pi(\sigma^d;\sigma^r)$ is constant in $\sigma^r$ there. -2. **$(c_t, z_t)$ representation**: emphasises that consumption is a martingale (random walk) - and that assets $b_t$ are encoded in consumption, so the impulse response function of - consumption is "box-shaped": a permanent shift in the level. +A meaningful comparison of rules must therefore move *off* the locus, for example by holding $\beta$ fixed at the benchmark and varying $\sigma^r$ alone. -We embedded this single-agent model in a Bewley equilibrium with a continuum of ex-post -heterogeneous consumers. +That experiment requires solving the full HST matrix problem rather than the scalar reduction of {ref}`rcs-scalar`, because once we leave the locus the marginal-utility process no longer summarises the decision rule. -The equilibrium gross interest rate $R = \beta^{-1}$ is supported by -constant average consumption, though the cross-section variance of consumption grows linearly with -age. +{doc}`robust_permanent_income` carries out that computation with the QuantEcon `LQ` and robust-control routines. -A complete-markets version of the same model achieves full risk sharing and a time-invariant -consumption distribution at the cost of more complex financial arrangements (Arrow securities). +## Concluding remarks -A concern for model misspecification, parameterised by $\sigma = -\theta^{-1} \leq 0$, alters the permanent income model. +A concern for model misspecification, parameterised by $\sigma = -\theta^{-1} \leq 0$, alters the permanent income model in ways that are subtle rather than dramatic. -A concern for robustness generates a precautionary savings motive even under quadratic preferences by distorting the conditional means of income shocks. +A concern for robustness generates a precautionary savings motive even under quadratic preferences, by distorting the conditional means of income shocks. -The distorted worst-case model makes the income process **more persistent**, shifting power toward low frequencies where the permanent income consumer is most vulnerable. +The distorted worst-case model makes the income process **more persistent**, shifting power toward the low frequencies where the permanent income consumer is most vulnerable. -The observational equivalence theorem {prf:ref}`thm-lqcs-oe1` shows that for quantities $(c_t, i_t)$ alone, a concern for robustness is indistinguishable from a reduction in $\beta$. +The observational equivalence theorem {prf:ref}`thm-rcs-oe1` shows that for quantities $(c_t, i_t)$ alone, a concern for robustness is indistinguishable from a reduction in $\beta$. -The reverse theorem {prf:ref}`thm-lqcs-oe2` shows that, starting from $\beta R = 1$, robustness is observationally equivalent to an *increase* in $\beta$, which imparts an upward drift to expected consumption. +The reverse theorem {prf:ref}`thm-rcs-oe2` shows that, starting from $\beta R = 1$, robustness is observationally equivalent to an *increase* in $\beta$, which imparts an upward drift to expected consumption. -Detection error probabilities provide a principled way to calibrate $\sigma$: choose $|\sigma|$ small enough that the approximating and worst-case models remain difficult to distinguish statistically. +Two disciplines bound how large a robustness concern can sensibly be. -The observationally equivalent $(\sigma, \hat\beta)$ pairs **do** have different implications for asset prices, a point explored further by HST in the asset-pricing context. +The breakdown point {eq}`eq:rcs-breakdown` is a hard mathematical limit beyond which no solution exists. -The robust Bewley economy shows how agents can have the same consumption decision rule and support the same equilibrium interest rate $R = \beta_0^{-1}$ while differing in their worst-case subjective income dynamics. +Detection error probabilities provide a softer and scale-free statistical discipline: choose $|\sigma|$ small enough that the approximating and worst-case models remain difficult to distinguish. + +The observationally equivalent $(\sigma, \hat\beta)$ pairs **do** have different implications for asset prices, a point pursued in {doc}`robust_permanent_income`. + +They also have different implications for the *beliefs* that agents hold, which is the subject of {doc}`lq_robust_bewley`. ## Exercises ```{exercise-start} -:label: lqcs_ex1 +:label: rcs_ex1 ``` -We translate from the benchmark Bewley economy to HST notation. - -Specialise the robust-control setup to the no-habit, no-capital LQ Bewley environment -($\lambda = \delta_h = 0$, $k_t = 0$), and let the endowment process be the two-factor model in -{eq}`eq:pi-twofactor`. +This exercise derives the breakdown point {eq}`eq:rcs-breakdown`. -1. Write the household state as $x_t = [a_t, z_t^\top]^\top$, where $a_t=-b_t$ is net assets, and derive matrices $(A, B, C)$ for the law of motion {eq}`eq:law0`. +1. Using {eq}`eq:rcs-zeta`, show that $\hat\beta(\sigma)\,\zeta(\sigma)^2 = \beta\,\zeta(\sigma)$. -2. Show that when $\sigma = 0$, the Bellman problem coincides with the LQ permanent-income - problem. +2. The minimizing agent's problem has a finite value only if $\hat\beta\zeta^2 < 1$. Use part 1 to show that this is equivalent to $\zeta < R$, and hence that + $\underline\sigma = -(1-\beta)^2/\alpha^2$. -3. Derive $\alpha^2$ and verify +3. Verify that $\hat\beta(\underline\sigma) = \beta^2$. -$$ -\alpha^2 = \sigma_1^2 + (1-\beta)^2\sigma_2^2. -$$ - -Interpret economically why the permanent and transitory components enter with different weights. +4. Explain why the breakdown point moves toward zero when the endowment becomes more volatile. ```{exercise-end} ``` -```{solution-start} lqcs_ex1 +```{solution-start} rcs_ex1 :class: dropdown ``` -Here is one solution: +Here is one solution. -1. With $x_t = [a_t, z_t^\top]^\top$ and budget law $a_{t+1} = R(a_t + y_t - c_t)$, $y_t = \check G z_t$, and $z_{t+1} = \check A z_t + \check C \epsilon_{t+1}$, the stacked law is +1. From {eq}`eq:rcs-zeta`, $\zeta = \beta/\hat\beta$, so $\hat\beta = \beta/\zeta$ and therefore $\hat\beta\zeta^2 = (\beta/\zeta)\zeta^2 = \beta\zeta$. -$$ -\begin{pmatrix} a_{t+1} \\ z_{t+1} \end{pmatrix} -= -\underbrace{\begin{pmatrix} R & R\check G \\ 0 & \check A \end{pmatrix}}_{A} -\begin{pmatrix} a_t \\ z_t \end{pmatrix} -+ -\underbrace{\begin{pmatrix} -R \\ 0 \end{pmatrix}}_{B} c_t -+ -\underbrace{\begin{pmatrix} 0 \\ \check C \end{pmatrix}}_{C}\epsilon_{t+1}. -$$ +2. By part 1 the condition $\hat\beta\zeta^2 < 1$ is $\beta\zeta < 1$, that is, $\zeta < \beta^{-1} = R$. - The sign of $B$ is negative because higher $c_t$ reduces asset accumulation $a_{t+1}$. + Substituting the closed form {eq}`eq:rcs-zeta` gives -2. At $\sigma=0$, the robust Bellman problem collapses to the ordinary LQ objective with no minimizing distortion term, so the planner/consumer problem is exactly the permanent-income problem with quadratic utility and linear constraints. + $$ + \left[1 + \frac{\sigma\alpha^2}{1-\beta}\right]^{-1} < \frac{1}{\beta} + \iff 1 + \frac{\sigma\alpha^2}{1-\beta} > \beta + \iff \sigma\alpha^2 > -(1-\beta)^2 , + $$ -3. From the $(c_t,z_t)$ representation, - -$$ -\Delta c_{t+1} = h\,\epsilon_{t+1}, -\qquad h = (1-\beta)\check G (I-\beta\check A)^{-1}\check C. -$$ + which is $\sigma > -(1-\beta)^2/\alpha^2 = \underline\sigma$. - In HST notation, $\alpha^2 = h h^\top$, and for the two-factor calibration $\check A=\mathrm{diag}(1,0)$ and $\check C=\mathrm{diag}(\sigma_1,\sigma_2)$, so +3. At $\sigma = \underline\sigma$ we have $\sigma\alpha^2 = -(1-\beta)^2$, so $\zeta = [1-(1-\beta)]^{-1} = \beta^{-1} = R$ and $\hat\beta = \beta/\zeta = \beta^2$. -$$ -\alpha^2 = \sigma_1^2 + (1-\beta)^2\sigma_2^2. -$$ +4. A more volatile endowment raises $\alpha^2$, and $\underline\sigma = -(1-\beta)^2/\alpha^2$ is therefore closer to zero. - Permanent shocks get unit weight because they shift lifetime resources one-for-one, while - transitory shocks are annuitised and therefore scaled by $(1-\beta)$ in consumption growth. + The economics is that $\sigma$ only ever matters through $\sigma\alpha^2$, so when the consumer faces more income risk a *smaller* $|\sigma|$ already buys a given amount of pessimism. ```{solution-end} ``` ```{exercise-start} -:label: lqcs_ex2 +:label: rcs_ex2 ``` -This exercise studies a continuum of robust but observationally equivalent Bewley consumers. - -Fix a benchmark pair $(\beta_0, \sigma = 0)$ with $R = \beta_0^{-1}$ and define +This exercise explains why a naive numerical solver for {eq}`eq:rcs-riccati` misbehaves. -$$ -\beta(\sigma) = \beta_0 + \frac{\sigma\alpha^2\beta_0}{1-\beta_0}, -\qquad \sigma \in [-\bar\sigma, 0]. -$$ - -Suppose a unit interval of consumers is indexed by $i$ with type $\sigma_i \in [-\bar\sigma, 0]$ -and discount factor $\beta_i = \beta(\sigma_i)$. +1. Verify algebraically that on the locus {eq}`eq:rcs-oe` the discriminant of {eq}`eq:rcs-riccati` equals $\bigl[(\sigma\alpha^2+(1-\beta)^2)/(1-\beta)\bigr]^2$. -1. Use {prf:ref}`thm-lqcs-oe1` to show that each type has the same consumption rule as the benchmark - $(\beta_0, 0)$ agent. +2. Conclude that the square root equals $|\sigma\alpha^2+(1-\beta)^2|/(1-\beta)$ and hence that the two roots are those in {eq}`eq:rcs-roots`. -2. Prove that aggregate consumption and bond-market clearing imply the same equilibrium interest - rate $R = \beta_0^{-1}$ as in the plain-vanilla Bewley model. +3. Write code that, for a grid of $\sigma$ in $(\underline\sigma,0)$, records which *sign* in front of the square root in {eq}`eq:rcs-riccati` yields the economically relevant root $P=-(1-\beta)^{-1}$. -3. Explain why agents can be observationally equivalent in quantities while still holding different - worst-case subjective models. + Confirm that the answer flips at $\sigma = \underline\sigma$. ```{exercise-end} ``` -```{solution-start} lqcs_ex2 +```{solution-start} rcs_ex2 :class: dropdown ``` -Here is one solution: +Here is one solution. -1. {prf:ref}`thm-lqcs-oe1` implies that if $(\sigma_i, \beta_i)$ lies on +1. With $u = \sigma\alpha^2$ and $\delta = 1-\beta$, equation {eq}`eq:rcs-oe` gives $\hat\beta - 1 = -\delta + u/\delta \cdot \beta/\beta$, more directly $\hat\beta-1+u = (u-\delta^2)/\delta$. -$$ -\beta_i = \beta_0 + \frac{\sigma_i\alpha^2\beta_0}{1-\beta_0}, -$$ + Hence the discriminant is - then type $i$ chooses the same decision rule as the benchmark $(0,\beta_0)$ agent and all types share the same consumption policy function $c_t = \mathcal C(a_t,z_t)$. + $$ + \frac{(u-\delta^2)^2}{\delta^2} + 4u + = \frac{(u-\delta^2)^2 + 4u\delta^2}{\delta^2} + = \frac{(u+\delta^2)^2}{\delta^2} . + $$ -2. Since all individual policy rules coincide with benchmark Bewley policies, aggregating over consumers gives the same goods- and bond-market clearing conditions and supports the same equilibrium $R=\beta_0^{-1}$. +2. Taking the square root gives $|u+\delta^2|/\delta$, and substituting the two signs into {eq}`eq:rcs-riccati` gives $P = -1/\delta$ and $P = \delta/u$. -3. Observational equivalence concerns quantities generated by optimal rules, so distinct $(\sigma_i,\beta_i)$ can generate the same $\{c_t^i,a_t^i\}$ while implying different internal worst-case beliefs. +3. The sign flips exactly where $u + \delta^2$ changes sign, that is at $u = -\delta^2$, which is $\sigma = \underline\sigma$. + +```{code-cell} ipython3 +σ_test = np.linspace(-0.01, 1.4 * σ_lo, 40) +signs = [] +for σ in σ_test: + roots, _ = solve_scalar_riccati(σ, β, α2) + signs.append('+' if np.argmin(np.abs(roots + 1 / (1 - β))) == 0 else '-') + +print(''.join(signs)) +flip = next(i for i in range(1, len(signs)) if signs[i] != signs[i - 1]) +print(f"sign flips between σ = {σ_test[flip - 1]:.5f} " + f"and σ = {σ_test[flip]:.5f}") +print(f"breakdown point σ̲ = {σ_lo:.5f}") +``` ```{solution-end} ``` ```{exercise-start} -:label: lqcs_ex3 +:label: rcs_ex3 ``` -This exercise separates quantities from beliefs without introducing an additional calibration. - -Consider two agents $a$ and $b$ in the robust Bewley economy with $\sigma^a < \sigma^b \leq 0$ and -$\beta^j = \beta_0 + \sigma^j\alpha^2\beta_0/(1-\beta_0)$ for $j \in \{a,b\}$. +This exercise calibrates $\sigma$ using the detection error probability. -1. Use {eq}`eq:bew_cinno` and {eq}`eq:bew_locus` to show that the two agents have the same consumption innovation $h\epsilon_{t+1}$. +Write a bisection that finds the $\sigma$ at which $\mathrm{DEP}(\sigma) = 0.2$, and report the associated $\hat\beta$, $\zeta$, and the ratio $\sigma/\underline\sigma$. -2. Show that if the two agents start from the same $(a_t,z_t)$ and observe the same shock - $\epsilon_{t+1}$, then their next-period choices of consumption and assets coincide. +Run it for $T = 40$ and $T = 160$. -3. Explain why the two agents can nevertheless disagree about the worst-case conditional mean of - $\epsilon_{t+1}$. - -Summarise what is and is not identified by data on quantities alone. +Take care: the target need not be attained on the admissible range $(\underline\sigma, 0]$, and your code should say so rather than silently return an endpoint. ```{exercise-end} ``` -```{solution-start} lqcs_ex3 +```{solution-start} rcs_ex3 :class: dropdown ``` -Here is one solution: +The DEP is decreasing in $|\sigma|$, so a bisection works, but we must first check that the target lies within reach. + +```{code-cell} ipython3 +def σ_for_target_dep(target, T, β, α2, tol=1e-5): + """ + Find σ ∈ (σ̲, 0) with DEP(σ) = target by bisection. -1. Equation {eq}`eq:bew_locus` places both agents on the observational-equivalence locus, so -{prf:ref}`thm-lqcs-oe1` implies that both use the benchmark consumption rule and therefore the same -innovation vector $h$ in {eq}`eq:bew_cinno`. + Returns None if the DEP never falls to the target on the admissible + range, which happens when the breakdown point binds before statistical + detectability does. + """ + α_loc = np.sqrt(α2) + lo, hi = 0.999 * (-(1 - β)**2 / α2), 0.0 # lo is the more negative end + + def dep_at(σ): + return detection_error_probability( + worst_case_persistence(σ, β, α2), α_loc, T=T) + + if dep_at(lo) > target: + return None + + while hi - lo > tol: + mid = 0.5 * (lo + hi) + if dep_at(mid) < target: + lo = mid # too easy to detect, move toward zero + else: + hi = mid + return 0.5 * (lo + hi) + + +for T in [40, 160]: + σ_star = σ_for_target_dep(0.2, T, β, α2) + if σ_star is None: + print(f"T = {T:>3}: DEP stays above 0.2 on the whole admissible " + f"range; the breakdown point binds first") + else: + ζ_star = worst_case_persistence(σ_star, β, α2) + print(f"T = {T:>3}: σ = {σ_star:.5f} β̂ = {β / ζ_star:.5f} " + f"ζ = {ζ_star:.5f} σ/σ̲ = {σ_star / σ_lo:.3f}") +``` -2. With a common state and common shock, both agents apply the same policy function and the same law -of motion, so $c_{t+1}^a=c_{t+1}^b$ and $a_{t+1}^a=a_{t+1}^b$. +A longer sample makes the two models easier to distinguish, so the $\sigma$ that keeps the DEP at $0.2$ moves closer to zero. -3. The minimizing feedback $K(\sigma^j,\beta^j)$ can differ across $j$, so the agents can attach -different worst-case conditional means to the same shock process even though their observable -choices coincide. +At $T = 40$ no admissible $\sigma$ has a DEP as low as $0.2$, so the breakdown point is the binding constraint. -Conclusion: quantities identify the equilibrium decision rule but not the decomposition between -impatience ($\beta$) and robustness ($\sigma$) along the observational-equivalence locus. +At $T = 160$ statistical detectability binds first, at about a quarter of the way to the breakdown point. ```{solution-end} ``` + +## Related lectures + +- {doc}`lq_permanent_income` develops the standard LQ permanent income model used here. +- {doc}`lq_bewley_complete_markets` studies the cross-section of consumption and the market structures that support it. +- {doc}`lq_robust_bewley` applies the observational-equivalence results of this lecture to build a Bewley economy with heterogeneous concerns about misspecification. +- {doc}`robust_permanent_income` treats risk-sensitive preferences, estimation, and asset pricing in the HST model. diff --git a/lectures/sargent_surico.md b/lectures/sargent_surico.md new file mode 100644 index 000000000..205d84bd1 --- /dev/null +++ b/lectures/sargent_surico.md @@ -0,0 +1,2096 @@ +--- +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.16.7 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +(sargent_surico)= +```{raw} jupyter +
+ + QuantEcon + +
+``` + +# Two Illustrations of the Quantity Theory of Money + +```{index} single: Quantity Theory of Money +``` + +```{contents} Contents +:depth: 2 +``` + +In addition to what's in Anaconda, this lecture uses `pandas_datareader` to download +macroeconomic data and `jax`, `numpyro` and `arviz` for the Hamiltonian Monte Carlo +section at the end: + +```{code-cell} ipython3 +:tags: [hide-output] + +!pip install pandas_datareader numpyro jax arviz +``` + +## Overview + +{cite:t}`Lucas1980` plotted long moving averages of U.S. inflation against long moving averages of U.S. money growth, and then long moving averages of a nominal interest rate against the same moving averages of money growth. + +Both scatter plots hugged a 45 degree line. + +Lucas read those two unit slopes as illustrating "two central implications of the quantity theory of money: that a given change in the rate of change in the quantity of money induces (i) an equal change in the rate of price inflation; and (ii) an equal change in nominal rates of interest." + +He also warned that a theory tells us "conditions under which one might expect them to break down." + +This lecture studies {cite:t}`SargentSurico2011`, which takes that warning seriously. + +The paper does three things. + +First, it extends Lucas's sample and shows that his two slopes are *not* stable across subperiods. + +Second, following {cite:t}`Whiteman1984`, it interprets a Lucas slope as the sum of coefficients in a two-sided distributed lag regression, an object that a time series model delivers as a ratio of spectral densities at frequency zero. + +Third, it estimates a small new Keynesian model on a pre-1984 sample and then perturbs *only* the monetary policy rule, showing that policy alone can move the two slopes across the range that the data display. + +We do all of this from scratch in Python. + +We write our own solver for linear rational expectations models, our own Kalman filter, and our own Metropolis-Hastings sampler, so that every step is visible. + +A final section then uses the estimated model as a test bed for Hamiltonian Monte Carlo, which turns out to require replacing the model solver with a differentiable one. + +Along the way we flag several places where the published paper's statement or implementation of its model needs care, and we check each of them numerically. + +Most are corrections, and one turns out to sharpen the paper's message rather than weaken it. + +Let's start with imports. + +```{code-cell} ipython3 +import datetime +import time +import warnings + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import pandas_datareader.data as web +from scipy import optimize, stats +from scipy.linalg import ordqz, solve_discrete_lyapunov, svd +from scipy.stats import invwishart + +warnings.filterwarnings('ignore') +plt.rcParams['figure.figsize'] = (10, 6) +plt.rcParams['axes.grid'] = True +plt.rcParams['grid.alpha'] = 0.3 +``` + +## Lucas's moving averages + +For a scalar series $x_t$ and $\beta \in [0, 1)$, {cite:t}`Lucas1980` formed the two-sided moving average + +$$ +\bar x_t(\beta) = a \sum_{k=-n}^{n} \beta^{|k|} x_{t+k}, +\qquad +a = \frac{(1-\beta)^2}{1 - \beta^2 - 2\beta^{n+1}(1-\beta)} , +$$ (eq:ss_filter) + +where $a$ makes the weights sum to one. + +Larger $\beta$ means a smoother series and a longer window. + +Setting $\beta = 0$ returns the raw data. + +We apply the weights to whatever observations are available, renormalizing so that they always sum to one. + +```{code-cell} ipython3 +def lucas_filter(x, beta=0.95): + """Two-sided exponentially weighted moving average of Lucas (1980).""" + x = np.asarray(x, dtype=float) + T = len(x) + out = np.empty(T) + for t in range(T): + w = beta ** np.abs(np.arange(T) - t) + out[t] = w @ x / w.sum() + return out +``` + +## Whiteman's reinterpretation + +{cite:t}`Whiteman1984` observed that fitting a straight line through a scatter of moving averages is an informal way of computing the *sum* of coefficients in a long two-sided distributed lag regression. + +Let $\{y_t, z_t\}$ be jointly covariance stationary with zero means, and project $y_t$ on the entire history and future of $z$, + +$$ +y_t = \sum_{j=-\infty}^{\infty} h_j z_{t-j} + \epsilon_t , +\qquad \mathbb{E}\, \epsilon_t z_{t-j} = 0 \ \ \forall j . +$$ (eq:ss_twosided) + +Write $S_y(\omega)$, $S_z(\omega)$ and $S_{yz}(\omega)$ for the spectral and cross-spectral densities and $h(\omega) = \sum_j h_j e^{-i\omega j}$ for the transfer function. + +Then + +$$ +h(\omega) = \frac{S_{yz}(\omega)}{S_z(\omega)} +\qquad \text{so} \qquad +\sum_{j=-\infty}^{\infty} h_j = h(0) = \frac{S_{yz}(0)}{S_z(0)} . +$$ (eq:ss_h0) + +Whiteman showed that as $\beta \to 1$ the regression coefficient of $\bar y_t(\beta)$ on $\bar z_t(\beta)$ converges to $h(0)$. + +So Lucas's two illustrations are claims that + +$$ +h_{\pi, \Delta m}(0) = 1 +\qquad \text{and} \qquad +h_{R, \Delta m}(0) = 1 , +$$ + +where $\pi$ is inflation, $\Delta m$ is money growth, and $R$ is a nominal interest rate. + +This is a useful reformulation because $h(0)$ is a *population* object that any time series model delivers. + +### From a state space model to $h(0)$ + +Suppose a model implies + +$$ +X_{t+1} = A X_t + B W_{t+1}, +\qquad +Y_{t+1} = C X_t + D W_{t+1}, +$$ (eq:ss_statespace) + +with $W_{t+1}$ a standard normal vector that is IID over time and $A$ stable. + +Iterating gives $X_t = \sum_{j \ge 0} A^j B W_{t-j}$, so the transfer function from $W$ to $Y$ is + +$$ +H(\zeta) = D + \zeta \, C (I - A\zeta)^{-1} B , +\qquad \zeta = e^{-i\omega}, +$$ + +and the spectral density matrix of $Y$ is $S_Y(\omega) = H(\zeta) H(\zeta)^{*} / (2\pi)$, that is + +$$ +2\pi S_Y(\omega) = +C (I - A e^{-i\omega})^{-1} B B' (I - A' e^{i\omega})^{-1} C' ++ D D' ++ e^{-i\omega} C (I - A e^{-i\omega})^{-1} B D' ++ e^{i\omega} D B' (I - A' e^{i\omega})^{-1} C' . +$$ (eq:ss_spectrum) + +```{note} +Equation (7) of {cite:t}`SargentSurico2011` reports only the first two terms of {eq}`eq:ss_spectrum`. + +The two cross terms drop out exactly when $B D' = 0$, that is when the shocks that move the state are orthogonal to the shocks that hit the observation equation. + +A one-line example shows that they matter otherwise: with scalars $A = 1/2$, $B = C = D = 1$, the correct $2\pi S_Y(0)$ is $(1 + 1/(1-1/2))^2 = 9$, while the truncated formula gives $1/(1-1/2)^2 + 1 = 5$. + +Nothing in what follows turns on this, because we write both our VAR and our solved model with $D = 0$, which is the natural way to write each of them. +``` + +At $\omega = 0$ with $D = 0$ everything collapses to a long-run multiplier. + +Writing $G = C(I - A)^{-1}B$, + +$$ +2\pi S_Y(0) = G G' , +\qquad +h_{y,z}(0) = \frac{[G G']_{yz}}{[G G']_{zz}} . +$$ (eq:ss_h0state) + +```{code-cell} ipython3 +def h_zero_from_state_space(A, B, C, iy, iz): + """h(0) for observable iy regressed on observable iz.""" + G = np.linalg.solve(np.eye(A.shape[0]) - A, B) + S0 = C @ G @ G.T @ C.T + return S0[iy, iz] / S0[iz, iz] +``` + +## Data + +{cite:t}`SargentSurico2011` splice FRED data onto the historical series of {cite:t}`BalkeGordon1986` to reach back to 1900. + +Those historical series are not on FRED, so we work with the post-1959 quarterly data that FRED does supply, and we extend the sample by twenty years beyond the paper's 2005 endpoint. + +That extension is the interesting part: it covers the balance sheet expansions after 2008 and after 2020, and the inflation of 2021-2023. + +We use M2 for money, the GDP deflator for prices, real GDP for output, and the three-month Treasury bill rate for the nominal interest rate. + +The paper's benchmark is a six-month commercial paper rate, which FRED no longer carries; its Appendix Table A reports that the three-month bill gives similar answers. + +Money growth, inflation, and output growth are 100 times quarterly log differences, and the interest rate is divided by four so that it too is a quarterly rate in percent. + +```{code-cell} ipython3 +start, end = datetime.datetime(1959, 1, 1), datetime.datetime(2030, 1, 1) +raw = web.DataReader(['M2SL', 'GDPDEF', 'GDPC1', 'TB3MS'], 'fred', start, end) + +q = pd.DataFrame({'M2': raw['M2SL'].resample('QS').mean(), + 'P': raw['GDPDEF'].resample('QS').mean(), + 'Y': raw['GDPC1'].resample('QS').mean(), + 'R': raw['TB3MS'].resample('QS').mean() / 4}).dropna() + +data = pd.DataFrame({'dm': 100 * np.log(q['M2']).diff(), + 'pi': 100 * np.log(q['P']).diff(), + 'dy': 100 * np.log(q['Y']).diff(), + 'R': q['R']}).dropna() + +VARS = ['dm', 'pi', 'R', 'dy'] +LABELS = {'dm': 'M2 growth', 'pi': 'inflation', + 'R': 'T-bill rate', 'dy': 'real GDP growth'} +print(f"sample: {data.index[0].date()} to {data.index[-1].date()}, T = {len(data)}") +data.describe().round(3) +``` + +Here are the raw series and their $\beta = 0.95$ moving averages. + +The shaded band marks 1959-1975, the closest we can come to Lucas's window. + +```{code-cell} ipython3 +filtered = pd.DataFrame({c: lucas_filter(data[c].values, 0.95) for c in VARS}, + index=data.index) + +fig, axes = plt.subplots(2, 2, figsize=(12, 7), sharex=True) +for ax, c in zip(axes.flat, VARS): + ax.plot(data.index, data[c], lw=0.8, color='0.65', label='raw') + ax.plot(filtered.index, filtered[c], lw=2.2, color='C0', + label=r'$\beta = 0.95$ filter') + ax.axvspan(pd.Timestamp('1959-01-01'), pd.Timestamp('1975-12-31'), + color='C1', alpha=0.12) + ax.set_title(LABELS[c]) + ax.legend(frameon=False, fontsize=9) +fig.suptitle('U.S. quarterly data and Lucas moving averages, percent per quarter') +plt.tight_layout() +plt.show() +``` + +Money growth, inflation, and the interest rate rise together into the early 1980s and then fall together. + +After 1984 that comovement is much harder to see, and the pandemic surge in M2 growth stands entirely apart from everything else. + +## Scatter plots + +Following the paper we plot second quarter observations of each year, which reduces the overlap between neighboring filtered values. + +```{code-cell} ipython3 +PERIODS = [('1959-1975', 1959, 1975), ('1960-1983', 1960, 1983), + ('1984-2007', 1984, 2007), ('2008-present', 2008, 2100)] + + +def scatter_panel(yvar, title): + q2 = filtered[filtered.index.quarter == 2] + pad = 0.15 * (q2[['dm', yvar]].values.max() - q2[['dm', yvar]].values.min()) + lim = np.array([q2[['dm', yvar]].values.min() - pad, + q2[['dm', yvar]].values.max() + pad]) + fig, axes = plt.subplots(1, 4, figsize=(14, 3.8), sharex=True, sharey=True) + for ax, (name, lo, hi) in zip(axes, PERIODS): + sub = q2[(q2.index.year >= lo) & (q2.index.year <= hi)] + x, y = sub['dm'].values, sub[yvar].values + ax.scatter(x, y, s=26, color='C0') + b = np.polyfit(x, y, 1) + xx = np.linspace(x.min(), x.max(), 2) + ax.plot(xx, np.polyval(b, xx), color='C3', lw=2, + label=f'slope = {b[0]:.2f}') + ax.plot(lim, lim, color='0.4', ls='--', lw=1, label='45 degrees') + ax.set(xlim=lim, ylim=lim, title=name, xlabel='money growth') + ax.legend(frameon=False, fontsize=8, loc='upper left') + axes[0].set_ylabel(title) + fig.suptitle(f'{title} against money growth, $\\beta = 0.95$ filtered') + plt.tight_layout() + plt.show() + + +scatter_panel('pi', 'inflation') +scatter_panel('R', 'interest rate') +``` + +Lucas's window and the pre-Volcker sample line up with the 45 degree line. + +The great moderation and the period since 2008 do not. + +## Regressions on filtered data + +Table 1 of the paper reports these slopes for a grid of $\beta$ values. + +```{code-cell} ipython3 +def filtered_slopes(betas=(0.95, 0.8, 0.5, 0.0)): + rows = {} + for beta in betas: + f = pd.DataFrame({c: lucas_filter(data[c].values, beta) for c in VARS}, + index=data.index) + for name, lo, hi in PERIODS + [('full sample', 1900, 2100)]: + m = ((f.index.year >= lo) & (f.index.year <= hi) + & (f.index.quarter == 2)) + x = f.loc[m, 'dm'].values + rows.setdefault(name, {})[('pi on dm', beta)] = np.polyfit(x, f.loc[m, 'pi'], 1)[0] + rows[name][('R on dm', beta)] = np.polyfit(x, f.loc[m, 'R'], 1)[0] + tab = pd.DataFrame(rows).T + order = [(v, b) for v in ['pi on dm', 'R on dm'] for b in betas] + return tab.reindex(columns=order) + + +filtered_slopes().round(2) +``` + +Two patterns stand out. + +The closer $\beta$ is to one, the larger the slopes, exactly as in Lucas's graphs and in the paper's Table 1. + +And the slopes fall steadily as we move through the sample, from above one in Lucas's window to negative after 2008. + +## Estimates of $h(0)$ from a VAR + +The filtered regressions are informal. + +Equation {eq}`eq:ss_h0state` lets us compute $h(0)$ properly from a fitted VAR. + +We fit a VAR(2) in money growth, inflation, the interest rate, and output growth, put it in companion form, and read $h(0)$ off the long-run multiplier matrix. + +We use a diffuse prior so that we can report posterior bands, as the paper does in its Figure 5. + +Draws whose companion matrix is explosive are discarded, since $h(0)$ is not defined for them. + +```{code-cell} ipython3 +def bvar_h0(Y, p=2, n_draws=500, seed=0): + """Posterior draws of h(0) from a VAR(p) under a diffuse prior.""" + rng = np.random.default_rng(seed) + T, n = Y.shape + X = np.column_stack([np.ones(T - p)] + [Y[p - l - 1:T - l - 1] for l in range(p)]) + Z, k = Y[p:], 1 + n * p + XXi = np.linalg.inv(X.T @ X) + B_hat = XXi @ X.T @ Z + E = Z - X @ B_hat + S, nu = E.T @ E, T - p - k + L_x = np.linalg.cholesky(XXi) + J = np.vstack([np.eye(n), np.zeros((n * (p - 1), n))]) + + out = np.full((n_draws, 2), np.nan) + for d in range(n_draws): + Sigma = invwishart.rvs(df=nu, scale=S, random_state=rng) + B = B_hat + L_x @ rng.standard_normal((k, n)) @ np.linalg.cholesky(Sigma).T + A = np.zeros((n * p, n * p)) + A[:n] = np.hstack([B[1 + l * n:1 + (l + 1) * n].T for l in range(p)]) + if p > 1: + A[n:, :n * (p - 1)] = np.eye(n * (p - 1)) + if np.max(abs(np.linalg.eigvals(A))) >= 1: + continue + G = np.linalg.solve(np.eye(n * p) - A, J) + S0 = G @ Sigma @ G.T + out[d] = [S0[1, 0] / S0[0, 0], S0[2, 0] / S0[0, 0]] + return out +``` + +```{code-cell} ipython3 +rows = [] +for name, lo, hi in PERIODS + [('full sample', 1900, 2100)]: + m = (data.index.year >= lo) & (data.index.year <= hi) + o = bvar_h0(data.loc[m, VARS].values) + qs = np.nanpercentile(o, [16, 50, 84], axis=0) + rows.append(dict(period=name, T=int(m.sum()), + h_pi=qs[1, 0], h_pi_lo=qs[0, 0], h_pi_hi=qs[2, 0], + h_R=qs[1, 1], h_R_lo=qs[0, 1], h_R_hi=qs[2, 1])) +pd.DataFrame(rows).set_index('period').round(2) +``` + +The subperiods are chosen by hand, which is the objection that {cite:t}`BoschenOtrok1994` raised against this style of evidence. + +So we also run the VAR through rolling twenty year windows. + +```{code-cell} ipython3 +years = np.arange(data.index.year.min(), data.index.year.max() - 18) +mid, lo68, hi68 = {}, {}, {} +for y in years: + m = (data.index.year >= y) & (data.index.year < y + 20) + if m.sum() < 70: + continue + o = bvar_h0(data.loc[m, VARS].values, n_draws=300, seed=int(y)) + qs = np.nanpercentile(o, [16, 50, 84], axis=0) + mid[y + 10], lo68[y + 10], hi68[y + 10] = qs[1], qs[0], qs[2] + +mid = pd.DataFrame(mid).T +lo68, hi68 = pd.DataFrame(lo68).T, pd.DataFrame(hi68).T + +fig, axes = plt.subplots(1, 2, figsize=(12, 4.2), sharex=True) +for j, (ax, ttl) in enumerate(zip(axes, [r'$h_{\pi,\Delta m}(0)$', + r'$h_{R,\Delta m}(0)$'])): + ax.plot(mid.index, mid[j], color='C3', lw=2, label='median') + ax.fill_between(mid.index, lo68[j], hi68[j], color='C3', alpha=0.18, + label='68% band') + ax.axhline(1, color='0.3', ls='--', lw=1, label='quantity theory') + ax.axhline(0, color='0.7', lw=0.8) + ax.set(title=ttl, xlabel='midpoint of 20 year window') + ax.legend(frameon=False, fontsize=9) +fig.suptitle('Low-frequency slopes from rolling VARs') +plt.tight_layout() +plt.show() +``` + +The two low-frequency slopes are high early in the sample, when Lucas wrote, and drift toward zero afterwards. + +One is rarely inside the 68 percent band after the mid 1980s. + +This is the instability that the rest of the lecture tries to explain. + +## A model for monetary policy analysis + +Section III of {cite:t}`SargentSurico2011` uses the log-linearized sticky price model of {cite:t}`Ireland2004`, which has price indexation, habit formation, and a unit root technology shock. + +The structure of the economy is + +$$ +\pi_t = \theta(1-\alpha_\pi) \mathbb{E}_t \pi_{t+1} + \theta \alpha_\pi \pi_{t-1} + + \kappa x_t - \tfrac{1}{\tau} e_t , +$$ (eq:ss_nkpc) + +$$ +x_t = (1-\alpha_x)\mathbb{E}_t x_{t+1} + \alpha_x x_{t-1} + - \sigma(R_t - \mathbb{E}_t \pi_{t+1}) + \sigma(1-\xi)(1-\rho_a) a_t , +$$ (eq:ss_is) + +$$ +\Delta m_t = \pi_t + z_t + \tfrac{1}{\sigma\gamma}\Delta x_t + - \tfrac{1}{\gamma}\Delta R_t + + \tfrac{1}{\gamma}(\Delta\chi_t - \Delta a_t) , +$$ (eq:ss_md) + +$$ +\tilde y_t = x_t + \xi a_t , +\qquad +\Delta y_t = \tilde y_t - \tilde y_{t-1} + z_t . +$$ (eq:ss_output) + +Here $\pi_t$ is inflation, $x_t$ the output gap, $\Delta m_t$ nominal money growth, $R_t$ the short rate, $\tilde y_t$ detrended output, and $z_t$ the growth rate of technology. + +Equation {eq}`eq:ss_nkpc` is a new Keynesian Phillips curve and {eq}`eq:ss_is` a new Keynesian IS curve. + +Equation {eq}`eq:ss_md` is the money demand relation of {cite:t}`McCallumNelson1999` and {cite:t}`Ireland2003`. + +The discount factor is $\theta$, $\alpha_\pi$ is indexation to past inflation, $\alpha_x$ is habit formation, $\kappa$ is the slope of the Phillips curve, $\sigma$ is the elasticity of intertemporal substitution, $\tau$ is the {cite:t}`Rotemberg1982` price adjustment cost, $\xi$ is the inverse Frisch elasticity, and $1/\gamma$ is the interest semi-elasticity of money demand. + +Four nonpolicy disturbances drive the economy: a markup shock $e_t$, a demand shock $a_t$, a money demand shock $\chi_t$, and technology $z_t$, + +$$ +e_t = \rho_e e_{t-1} + \varepsilon_{et}, +\quad +a_t = \rho_a a_{t-1} + \varepsilon_{at}, +\quad +\chi_t = \rho_\chi \chi_{t-1} + \varepsilon_{\chi t}, +\quad +z_t = \varepsilon_{zt} . +$$ (eq:ss_shocks) + +Monetary policy is either a money growth rule + +$$ +\Delta m_t = \rho_m \Delta m_{t-1} + (1-\rho_m)(\phi_\pi \pi_t + \phi_x x_t) + + \varepsilon_{mt} +$$ (eq:ss_mrule) + +or a Taylor rule + +$$ +R_t = \rho_r R_{t-1} + (1-\rho_r)(\psi_\pi \pi_t + \psi_x x_t) + \varepsilon_{Rt} . +$$ (eq:ss_trule) + +The observables are $[\Delta m_t, \pi_t, R_t, \Delta y_t]$. + +## Three things worth checking + +Before estimating anything it pays to look at what this parameterization can and cannot deliver. + +### $\tau$ is not identified + +The price adjustment cost $\tau$ enters the model in exactly one place, the term $e_t/\tau$ in {eq}`eq:ss_nkpc`. + +Since $e_t$ is an AR(1) with a freely estimated innovation standard deviation $\sigma_e$, the process $e_t/\tau$ is an AR(1) with innovation standard deviation $\sigma_e/\tau$. + +So the likelihood depends on $(\tau, \sigma_e)$ only through the ratio $\sigma_e/\tau$. + +We verify this numerically below. + +Table 2 of the paper nevertheless reports a posterior for $\tau$ with mean 3.51 and a 5-95 interval of $[1.99, 4.99]$, against a prior with mean 4 and interval $[2.51, 5.77]$. + +That posterior is not evidence about the Rotemberg adjustment cost. + +Its leftward shift is what the prior on $\sigma_e$ implies once the data pin down the ratio: the data want $\sigma_e/\tau \approx 0.31$, while the inverse gamma prior on $\sigma_e$ has mean $0.3$ and pulls $\sigma_e$ down, dragging $\tau$ with it. + +Because $\tau$ contributes nothing, we fix it at its prior mean and estimate one fewer parameter. + +Nothing in the fit, or in $h(0)$, changes. + +### The first unit slope is an identity when money is exogenous + +Group the terms of the money demand equation {eq}`eq:ss_md` as + +$$ +\Delta m_t = \pi_t + z_t + \Delta v_t , +\qquad +v_t = \tfrac{1}{\sigma\gamma} x_t - \tfrac{1}{\gamma} R_t + + \tfrac{1}{\gamma}(\chi_t - a_t) . +$$ (eq:ss_qtm) + +Everything other than $\pi_t + z_t$ is a *first difference*. + +The filter $1 - e^{-i\omega}$ vanishes at $\omega = 0$, so $\Delta v_t$ contributes nothing to any spectral density or cross spectrum at frequency zero. + +Therefore + +$$ +h_{\pi, \Delta m}(0) += \frac{S_{\pi,\pi+z}(0)}{S_{\pi+z}(0)} += 1 - \frac{S_{z,\pi+z}(0)}{S_{\pi+z}(0)} . +$$ (eq:ss_decomp) + +This yields a proposition. + +```{prf:proposition} +:label: ss_prop + +Under the money growth rule {eq}`eq:ss_mrule` with $\phi_\pi = \phi_x = 0$, + +$$ +h_{\pi, \Delta m}(0) = 1 +$$ + +exactly, whatever the values of the other parameters. +``` + +```{prf:proof} +With $\phi_\pi = \phi_x = 0$ the rule reads $\Delta m_t = \rho_m \Delta m_{t-1} + \varepsilon_{mt}$, so money growth is driven by $\varepsilon_m$ alone. + +From {eq}`eq:ss_qtm`, $\pi_t + z_t = \Delta m_t - \Delta v_t$, and $\Delta v$ contributes nothing at frequency zero, so $S_{z,\pi+z}(0) = S_{z,\Delta m}(0)$. + +Technology is exogenous and orthogonal to $\varepsilon_m$, so $S_{z, \Delta m}(\omega) \equiv 0$. + +Now apply {eq}`eq:ss_decomp`. +``` + +This is worth thinking about. + +Lucas's *first* illustration prevails whenever money growth is econometrically exogenous, which is the case {cite:t}`Whiteman1984` assumed. + +It is a property of the money demand specification, not of a deep monetary neutrality. + +Departures from one arise only through $S_{z,\pi+z}(0)$, and that term is nonzero only when the policy rule makes money growth respond to the endogenous variables, so that technology becomes correlated with money growth at frequency zero. + +The *second* illustration has no such backing. + +The same argument gives $h_{R,\Delta m}(0) = S_{R,\pi+z}(0)/S_{\pi+z}(0)$, which equals one only if the low-frequency behavior of the ex ante real rate happens to cooperate. + +So within this structural model, Lucas's two illustrations are not on the same footing. + +### The reported posterior for $\alpha_x$ stops at one half + +Table 2 reports a posterior for the habit parameter $\alpha_x$ with mean 0.4775 and a 95th percentile of exactly 0.5000, under a beta prior on $[0,1]$. + +An interval that ends on a round number is usually a constraint rather than a finding. + +It is not a determinacy constraint: we check below that the model has a unique stable solution on both sides of $\alpha_x = 0.5$. + +We therefore let $\alpha_x$ range over the whole unit interval. + +## Solving the model + +The model is a linear rational expectations system, which we solve with the method of {cite:t}`Sims2002gensys`. + +Write it in the canonical form + +$$ +\Gamma_0 y_t = \Gamma_1 y_{t-1} + \Psi \varepsilon_t + \Pi \eta_t , +$$ (eq:ss_gensys) + +where $\varepsilon_t$ collects the structural shocks and $\eta_t$ collects expectational errors, one for each variable whose expectation appears in the system. + +The trick is to add $\mathbb{E}_t \pi_{t+1}$ and $\mathbb{E}_t x_{t+1}$ to the vector of variables, together with the identities + +$$ +\pi_t = \mathbb{E}_{t-1}\pi_t + \eta^\pi_t , +\qquad +x_t = \mathbb{E}_{t-1}x_t + \eta^x_t . +$$ + +The solution algorithm takes a generalized Schur decomposition of the matrix pencil $(\Gamma_0, \Gamma_1)$, orders the stable generalized eigenvalues first, and then asks whether the expectational errors can be chosen to kill the explosive block. + +If they can, a stable solution exists; if they can be chosen in only one way, it is unique. + +```{code-cell} ipython3 +SMALL = 1e-6 + + +def gensys(g0, g1, psi, pi, div=1.01): + """ + Solve g0 @ y[t] = g1 @ y[t-1] + psi @ eps[t] + pi @ eta[t]. + + Returns (G1, impact, eu) so that, when eu == (1, 1), + + y[t] = G1 @ y[t-1] + impact @ eps[t] + + is the unique stable solution. eu[0] flags existence, eu[1] uniqueness. + """ + n = g0.shape[0] + S, T, alpha, beta, Q, Z = ordqz(g0, g1, output='complex', + sort=lambda a, b: abs(b) < div * abs(a)) + nunstab = int(np.sum(abs(beta) >= div * abs(alpha))) + ns = n - nunstab + if np.any((abs(alpha) < SMALL) & (abs(beta) < SMALL)): + return None, None, (-2, -2) + + q = Q.conj().T + q1, q2 = q[:ns], q[ns:] + + def trimmed_svd(M): + u, d, vh = svd(M) + keep = d > SMALL + r = len(d) + return u[:, :r][:, keep], d[keep], vh.conj().T[:, :r][:, keep] + + u2, d2, v2 = trimmed_svd(q2 @ pi) # can eta kill the explosive block? + u1, d1, v1 = trimmed_svd(q1 @ pi) # is the choice unique? + + exist = len(d2) >= nunstab + if v1.shape[1] == 0: + unique = True + else: + loose = v1 - v2 @ v2.conj().T @ v1 + unique = np.sum(svd(loose, compute_uv=False) > SMALL * n) == 0 + eu = (int(exist), int(unique)) + if not exist: + return None, None, eu + + W = u2 @ np.diag(1 / d2) @ v2.conj().T @ v1 @ np.diag(d1) @ u1.conj().T + tmat = np.hstack([np.eye(ns), -W.conj().T]) + G0 = np.vstack([tmat @ S, + np.hstack([np.zeros((nunstab, ns)), np.eye(nunstab)])]) + G0i = np.linalg.inv(G0) + G1 = np.real(Z @ (G0i @ np.vstack([tmat @ T, np.zeros((nunstab, n))])) + @ Z.conj().T) + impact = np.real(Z @ (G0i @ np.vstack([tmat @ q @ psi, + np.zeros((nunstab, psi.shape[1]))]))) + return G1, impact, eu +``` + +Before trusting it, we check it against a model whose solution we can write down. + +In the textbook three-equation new Keynesian model + +$$ +\pi_t = \theta \mathbb{E}_t \pi_{t+1} + \kappa x_t, +\quad +x_t = \mathbb{E}_t x_{t+1} - \sigma(R_t - \mathbb{E}_t \pi_{t+1}) + g_t, +\quad +R_t = \phi_\pi \pi_t , +$$ + +with $g_t$ an AR(1), the equilibrium is $\pi_t = \pi_g g_t$ and $x_t = x_g g_t$ where $(\pi_g, x_g)$ solves a two by two linear system, and the equilibrium is unique if and only if $\phi_\pi > 1$. + +```{code-cell} ipython3 +def toy_nk(phi_pi, theta=0.99, kappa=0.1, sigma=1.0, rho=0.8): + """y = [pi, x, R, g, E pi(+1), E x(+1)].""" + g0, g1 = np.zeros((6, 6)), np.zeros((6, 6)) + psi, pie = np.zeros((6, 1)), np.zeros((6, 2)) + g0[0, 0], g0[0, 4], g0[0, 1] = 1, -theta, -kappa + g0[1, 1], g0[1, 5], g0[1, 2], g0[1, 4], g0[1, 3] = 1, -1, sigma, -sigma, -1 + g0[2, 2], g0[2, 0] = 1, -phi_pi + g0[3, 3], g1[3, 3], psi[3, 0] = 1, rho, 1 + g0[4, 0], g1[4, 4], pie[4, 0] = 1, 1, 1 + g0[5, 1], g1[5, 5], pie[5, 1] = 1, 1, 1 + return g0, g1, psi, pie + + +for phi in [0.5, 1.5, 3.0]: + G1, impact, eu = gensys(*toy_nk(phi)) + msg = 'unique stable solution' if eu == (1, 1) else 'indeterminate' + print(f'phi_pi = {phi}: eu = {eu} ({msg})') + +theta, kappa, sigma, rho, phi = 0.99, 0.1, 1.0, 0.8, 1.5 +M = np.array([[1 - theta * rho, -kappa], [sigma * (phi - rho), 1 - rho]]) +exact = np.linalg.solve(M, np.array([0.0, 1.0])) +G1, impact, _ = gensys(*toy_nk(phi)) +print(f'\nanalytic (pi_g, x_g) = {np.round(exact, 8)}') +print(f'gensys (pi_g, x_g) = {np.round(impact[:2, 0], 8)}') +``` + +The solver reproduces the analytic solution to machine precision and correctly refuses to deliver one when $\phi_\pi < 1$. + +## The Sargent-Surico model in canonical form + +Now we write {eq}`eq:ss_nkpc` through {eq}`eq:ss_trule` in the form {eq}`eq:ss_gensys`. + +The vector of variables is + +$$ +y_t = [\pi_t,\ x_t,\ \Delta m_t,\ R_t,\ e_t,\ a_t,\ \chi_t,\ z_t,\ + \mathbb{E}_t \pi_{t+1},\ \mathbb{E}_t x_{t+1}]' , +$$ + +and the shocks are $\varepsilon_t = [\varepsilon_{et}, \varepsilon_{at}, \varepsilon_{\chi t}, \varepsilon_{zt}, \varepsilon_{mt}]'$. + +```{code-cell} ipython3 +PI, X, DM, R, E, A_, CH, Z, EPI, EX = range(10) +NY = 10 + + +def canonical(p, rule='money'): + """Matrices g0, g1, psi, pi of equations (SS-NKPC) through (SS-Taylor).""" + th, api, kap, tau = p['theta'], p['alpha_pi'], p['kappa'], p['tau'] + ax, sig, xi, gam = p['alpha_x'], p['sigma'], p['xi'], p['gamma'] + g0, g1 = np.zeros((NY, NY)), np.zeros((NY, NY)) + psi, pie = np.zeros((NY, 5)), np.zeros((NY, 2)) + + # Phillips curve + g0[0, PI], g0[0, EPI], g0[0, X], g0[0, E] = 1, -th * (1 - api), -kap, 1 / tau + g1[0, PI] = th * api + # IS curve + g0[1, X], g0[1, EX], g0[1, R], g0[1, EPI] = 1, -(1 - ax), sig, -sig + g0[1, A_] = -sig * (1 - xi) * (1 - p['rho_a']) + g1[1, X] = ax + # money demand + g0[2, DM], g0[2, PI], g0[2, Z] = 1, -1, -1 + for col, coef in [(X, -1 / (sig * gam)), (R, 1 / gam), + (CH, -1 / gam), (A_, 1 / gam)]: + g0[2, col], g1[2, col] = coef, coef + # policy rule + if rule == 'money': + rm = p['rho_m'] + g0[3, DM] = 1 + g0[3, PI], g0[3, X] = -(1 - rm) * p['phi_pi'], -(1 - rm) * p['phi_x'] + g1[3, DM] = rm + else: + rr = p['rho_r'] + g0[3, R] = 1 + g0[3, PI], g0[3, X] = -(1 - rr) * p['psi_pi'], -(1 - rr) * p['psi_x'] + g1[3, R] = rr + psi[3, 4] = 1 + # exogenous shocks + for row, (v, rho, k) in enumerate([(E, p['rho_e'], 0), (A_, p['rho_a'], 1), + (CH, p['rho_chi'], 2), (Z, 0.0, 3)], 4): + g0[row, v], g1[row, v], psi[row, k] = 1, rho, 1 + # expectational identities + g0[8, PI], g1[8, EPI], pie[8, 0] = 1, 1, 1 + g0[9, X], g1[9, EX], pie[9, 1] = 1, 1, 1 + return g0, g1, psi, pie +``` + +Output growth in {eq}`eq:ss_output` needs $x_{t-1}$ and $a_{t-1}$, so the state carries those two lags, + +$$ +S_t = [y_t',\ x_{t-1},\ a_{t-1}]' , +\qquad +S_t = A S_{t-1} + B \varepsilon_t , +\qquad +Y_t = C S_t . +$$ + +There is no measurement error, so $D = 0$ and {eq}`eq:ss_h0state` applies directly. + +The following issue arises here. + +The `div` argument of `gensys` splits stable from explosive generalized eigenvalues, and it has to sit a little above one for the split to be numerically reliable. + +So a solution that `gensys` reports as unique and stable can still have a root of the transition matrix at or just above one. + +Both $h(0)$ and the Kalman filter need covariance stationarity, so we test for it explicitly rather than trusting the flags. + +```{code-cell} ipython3 +def state_space(p, rule='money'): + """ + Return A, B, C for S[t] = A S[t-1] + B eps[t] and Y[t] = C S[t], + together with a status string that is 'ok' when the equilibrium is + unique and covariance stationary. + """ + G1, impact, eu = gensys(*canonical(p, rule)) + if eu[0] != 1: + return None, None, None, 'no stable solution' + if eu[1] != 1: + return None, None, None, 'indeterminate' + sd = np.diag([p['sig_e'], p['sig_a'], p['sig_chi'], p['sig_z'], p['sig_m']]) + n = NY + 2 + A, B = np.zeros((n, n)), np.zeros((n, 5)) + A[:NY, :NY] = G1 + A[NY, X], A[NY + 1, A_] = 1, 1 + B[:NY] = impact @ sd + C = np.zeros((4, n)) # [dm, pi, R, dy] + C[0, DM], C[1, PI], C[2, R] = 1, 1, 1 + C[3, X], C[3, NY] = 1, -1 + C[3, A_], C[3, NY + 1] = p['xi'], -p['xi'] + C[3, Z] = 1 + if np.max(np.abs(np.linalg.eigvals(A))) > 1 - 1e-9: + return None, None, None, 'unit or explosive root' + return A, B, C, 'ok' + + +def h_zero(p, rule='money'): + """(h_pi,dm(0), h_R,dm(0)) implied by the model at parameters p.""" + A, B, C, status = state_space(p, rule) + if status != 'ok': + return np.nan, np.nan + return (h_zero_from_state_space(A, B, C, 1, 0), + h_zero_from_state_space(A, B, C, 2, 0)) +``` + +This test matters. + +Without it, a point like $\phi_\pi = 1$, $\phi_x = 0$ returns $h_{\pi,\Delta m}(0) = h_{R,\Delta m}(0) = 1.000$, an apparently perfect confirmation of both of Lucas's illustrations. + +The transition matrix there has a root of exactly one, the spectral density at frequency zero does not exist, and those two numbers come from inverting a matrix whose condition number is about $10^{16}$. + +That point sits in the upper right corner of the range that Figure 6 of the paper plots. + + +We show below that at the paper's own Table 2 posterior means a whole strip of that range has no covariance stationary equilibrium. + +We can now check our implementation against the paper. + +At the posterior means of Table 2, {cite:t}`SargentSurico2011` report implied values $h_{\pi,\Delta m}(0) = 1.0068$ and $h_{R,\Delta m}(0) = 0.8163$. + +```{code-cell} ipython3 +PAPER = dict(theta=0.9901, alpha_pi=0.8815, kappa=0.0324, tau=3.5126, + alpha_x=0.4775, sigma=0.0997, xi=3.0319, gamma=3.8128, + phi_pi=0.2312, phi_x=-0.1971, rho_m=0.7428, rho_e=0.5645, + rho_a=0.9241, rho_chi=0.5024, sig_e=1.0922, sig_a=0.7226, + sig_chi=0.2388, sig_z=1.5845, sig_m=1.1457) + +print('h at the posterior means of Table 2: %.4f, %.4f' % h_zero(PAPER)) +print('reported in Table 2: 1.0068, 0.8163') +``` + +Our independent implementation lands on the paper's numbers. + +That gives us confidence that we have read equations {eq}`eq:ss_nkpc` through {eq}`eq:ss_mrule` the way their authors intended. + +### Checking the three claims + +Now we verify the three observations of the previous section. + +```{code-cell} ipython3 +print('scaling tau and sigma_e together leaves h(0) untouched:') +for f in [0.5, 1.0, 2.0, 8.0]: + q = dict(PAPER, tau=PAPER['tau'] * f, sig_e=PAPER['sig_e'] * f) + print(f' tau, sigma_e scaled by {f:4.1f}: h_pi = {h_zero(q)[0]:.10f}') + +print('\nwith phi_pi = phi_x = 0, h_pi is exactly one for any other parameters:') +rng = np.random.default_rng(0) +for i in range(5): + q = dict(PAPER, phi_pi=0.0, phi_x=0.0, + kappa=rng.uniform(0.01, 0.3), alpha_pi=rng.uniform(0.1, 0.9), + alpha_x=rng.uniform(0.1, 0.9), sigma=rng.uniform(0.05, 0.5), + gamma=rng.uniform(1, 8), rho_m=rng.uniform(0, 0.95), + sig_z=rng.uniform(0.2, 3.0), xi=rng.uniform(0.5, 5.0)) + print(f' draw {i}: h_pi = {h_zero(q)[0]:.12f} h_R = {h_zero(q)[1]:.4f}') + +print('\ndeterminacy on both sides of alpha_x = 0.5:') +for ax in [0.30, 0.45, 0.499, 0.501, 0.60, 0.80]: + print(f' alpha_x = {ax:5.3f}: {state_space(dict(PAPER, alpha_x=ax))[3]}') +``` + +All three claims hold. + +The likelihood is invariant to scaling $(\tau, \sigma_e)$, the first unit slope is an identity under exogenous money growth, and the model is determinate on both sides of $\alpha_x = 1/2$. + +Notice also that $h_{R,\Delta m}(0)$ moves all over the place across those same draws, which is the asymmetry between the two illustrations that {prf:ref}`ss_prop` predicts. + +Finally, here is the strip of the Figure 6 policy grid on which, at the paper's own posterior means, no covariance stationary equilibrium exists. + +```{code-cell} ipython3 +for phi_x in [0.0, -0.25, -0.5]: + bad = [round(g, 2) for g in np.linspace(-3, 1, 41) + if state_space(dict(PAPER, phi_pi=g, phi_x=phi_x))[3] != 'ok'] + print(f'phi_x = {phi_x:5.2f}: no stationary equilibrium at phi_pi in ' + f'{bad if bad else "none of the grid"}') +``` + +The strip is far from the estimated policy rule, so the paper's conclusions are not at stake. + +But an $h(0)$ reported over it would be a number computed from a spectral density that does not exist. + +## Bayesian estimation + +We estimate the money growth rule version on 1960:I-1983:IV, the paper's sample. + +The four observables are demeaned, since the model is written in deviations from a steady state. + +```{code-cell} ipython3 +mask = (data.index.year >= 1960) & (data.index.year <= 1983) +Y_est = data.loc[mask, VARS].values +Y_est = Y_est - Y_est.mean(0) +print(f'estimation sample: {mask.sum()} quarters, ' + f'{data.index[mask][0].date()} to {data.index[mask][-1].date()}') +``` + +### The likelihood + +Given $(A, B, C)$ the likelihood follows from the Kalman filter. + +We initialize the state at its unconditional mean and covariance, the latter solving the discrete Lyapunov equation $P = A P A' + BB'$. + +```{code-cell} ipython3 +def loglik(p, Y, rule='money'): + """Kalman filter log likelihood of Y (T x 4) at parameters p.""" + A, B, C, status = state_space(p, rule) + if status != 'ok': + return -np.inf + Q = B @ B.T + try: + P = solve_discrete_lyapunov(A, Q) + except Exception: + return -np.inf + s = np.zeros(A.shape[0]) + ll, const = 0.0, Y.shape[1] * np.log(2 * np.pi) + for t in range(Y.shape[0]): + s = A @ s + P = A @ P @ A.T + Q + v = Y[t] - C @ s # forecast error + PCt = P @ C.T + F = C @ PCt # its covariance + try: + L = np.linalg.cholesky(F) + except np.linalg.LinAlgError: + return -np.inf + u = np.linalg.solve(L, v) + ll -= 0.5 * (const + 2 * np.sum(np.log(np.diag(L))) + u @ u) + K = np.linalg.solve(F, PCt.T).T # Kalman gain + s, P = s + K @ v, P - K @ PCt.T + return ll if np.isfinite(ll) else -np.inf +``` + +### Priors + +We use the prior means and standard deviations of Table 2, matching moments to pick the parameters of each family. + +```{code-cell} ipython3 +def beta_prior(m, s): + nu = m * (1 - m) / s ** 2 - 1 + return stats.beta(m * nu, (1 - m) * nu) + + +def gamma_prior(m, s): + return stats.gamma(m ** 2 / s ** 2, scale=s ** 2 / m) + + +def invgamma_prior(m, s): + a = m ** 2 / s ** 2 + 2 + return stats.invgamma(a, scale=m * (a - 1)) + + +PRIOR, SUPPORT = {}, {} +for n, (m, s) in [('theta', (0.99, 0.005)), ('alpha_pi', (0.5, 0.2)), + ('alpha_x', (0.5, 0.2)), ('rho_m', (0.5, 0.05)), + ('rho_e', (0.5, 0.1)), ('rho_a', (0.5, 0.1)), + ('rho_chi', (0.5, 0.1))]: + PRIOR[n], SUPPORT[n] = beta_prior(m, s), (1e-6, 1 - 1e-6) +for n, (m, s) in [('kappa', (0.3, 0.1)), ('sigma', (0.1, 0.05)), + ('xi', (2.0, 1.0)), ('gamma', (4.0, 1.0))]: + PRIOR[n], SUPPORT[n] = gamma_prior(m, s), (1e-8, np.inf) +for n in ['phi_pi', 'phi_x']: + PRIOR[n], SUPPORT[n] = stats.norm(0, 0.5), (-np.inf, np.inf) +for n in ['sig_e', 'sig_a', 'sig_chi', 'sig_z', 'sig_m']: + PRIOR[n], SUPPORT[n] = invgamma_prior(0.3, 1.0), (1e-8, np.inf) + +FREE = ['theta', 'alpha_pi', 'kappa', 'alpha_x', 'sigma', 'xi', 'gamma', + 'phi_pi', 'phi_x', 'rho_m', 'rho_e', 'rho_a', 'rho_chi', + 'sig_e', 'sig_a', 'sig_chi', 'sig_z', 'sig_m'] +TAU = 4.0 # not identified; fixed at its prior mean + + +def unpack(v): + return dict({n: float(x) for n, x in zip(FREE, v)}, tau=TAU) + + +def log_post(v, Y): + lp = 0.0 + for n, x in zip(FREE, v): + lo, hi = SUPPORT[n] + if not lo < x < hi: + return -np.inf + lp += PRIOR[n].logpdf(x) + return lp + loglik(unpack(v), Y) +``` + +### The posterior mode + +We start the search at the paper's posterior means, use Powell's derivative free method, and polish with a gradient method. + +```{code-cell} ipython3 +v_start = np.array([PAPER[n] for n in FREE]) +neg_log_post = lambda v: -log_post(v, Y_est) + +res = optimize.minimize(neg_log_post, v_start, method='Powell', + options=dict(maxiter=20000, maxfev=20000)) +res = optimize.minimize(neg_log_post, res.x, method='L-BFGS-B', + bounds=[SUPPORT[n] for n in FREE]) +v_mode = res.x +print(f'log posterior at the paper\'s means: {log_post(v_start, Y_est):10.3f}') +print(f'log posterior at the mode: {-res.fun:10.3f}') +print('h(0) at the mode: %.4f, %.4f' % h_zero(unpack(v_mode))) +``` + +```{note} +Starting the search at the paper's estimates finds the mode nearest to them. + +The Hamiltonian Monte Carlo section at the end of this lecture discovers that this is +a *local* mode, and that the posterior has a second one with higher density. +``` + +### Random walk Metropolis-Hastings + +The standard proposal covariance is the inverse Hessian of the negative log posterior at the mode, which we compute by finite differences. + +```{code-cell} ipython3 +def numerical_hessian(f, v, rel=1e-4): + n = len(v) + h = rel * np.maximum(np.abs(v), 1e-2) + H = np.zeros((n, n)) + for i in range(n): + for j in range(i, n): + ei, ej = np.zeros(n), np.zeros(n) + ei[i], ej[j] = h[i], h[j] + H[i, j] = H[j, i] = (f(v + ei + ej) - f(v + ei - ej) + - f(v - ei + ej) + f(v - ei - ej)) / (4 * h[i] * h[j]) + return H + + +H = numerical_hessian(neg_log_post, v_mode) +Sigma_prop = np.linalg.inv(H) +print('proposal covariance is positive definite:', + np.all(np.linalg.eigvalsh(H) > 0)) +``` + +```{code-cell} ipython3 +def rwmh(Y, v0, Sigma, n_draws, c=0.45, seed=42): + """Random walk Metropolis-Hastings in the natural parameter space.""" + rng = np.random.default_rng(seed) + L = np.linalg.cholesky(Sigma) + v, lp = v0.copy(), log_post(v0, Y) + draws, n_acc = np.empty((n_draws, len(v))), 0 + for i in range(n_draws): + cand = v + c * (L @ rng.standard_normal(len(v))) + lp_cand = log_post(cand, Y) + if np.log(rng.random()) < lp_cand - lp: + v, lp, n_acc = cand, lp_cand, n_acc + 1 + draws[i] = v + return draws, n_acc / n_draws +``` + +Draws that fall outside the support of a prior get $-\infty$ and are rejected, so we can work in the natural parameter space and skip transformations. + +The chain below is short enough to run while you read; a serious application would use many more draws. + +```{code-cell} ipython3 +N_DRAWS, BURN = 30_000, 10_000 + +t0 = time.time() +draws, acc_rate = rwmh(Y_est, v_mode, Sigma_prop, N_DRAWS) +rwmh_seconds = time.time() - t0 + +kept = draws[BURN::5] +print(f'acceptance rate {acc_rate:.3f}, {len(kept)} retained draws, ' + f'{rwmh_seconds:.0f} seconds') +``` + +```{code-cell} ipython3 +fig, axes = plt.subplots(2, 3, figsize=(12, 5)) +for ax, n in zip(axes.flat, ['phi_pi', 'phi_x', 'rho_m', + 'alpha_pi', 'kappa', 'sigma']): + j = FREE.index(n) + ax.plot(draws[:, j], lw=0.4, color='C0') + ax.axvline(BURN, color='C3', ls='--', lw=1) + ax.set_title(n) +fig.suptitle('Metropolis-Hastings traces, dashed line ends the burn-in') +plt.tight_layout() +plt.show() +``` + +```{code-cell} ipython3 +summary = pd.DataFrame({ + 'prior mean': [PRIOR[n].mean() for n in FREE], + 'post. mean': kept.mean(0), + '5th': np.percentile(kept, 5, axis=0), + '95th': np.percentile(kept, 95, axis=0), + 'paper': [PAPER[n] for n in FREE]}, index=FREE) +summary.round(4) +``` + +Our estimates are recognizably those of Table 2, with the differences one expects from a different interest rate series and different data vintages. + +The Phillips curve is strongly backward looking and flat, the IS curve much less so. + +Most important for what follows, the money growth rule responds weakly to inflation, with a posterior for $\phi_\pi$ that straddles zero, and it carries substantial smoothing. + +That is the paper's central reading of the pre-1984 regime: the Federal Reserve put persistent, nearly exogenous, movement into money growth. + +By {prf:ref}`ss_prop` that is precisely the configuration in which Lucas's first illustration holds, which is also why the paper's own posterior for $h_{\pi,\Delta m}(0)$ in Table 2 is so remarkably tight around one. + +### The posterior distribution of $h(0)$ + +Every draw of the structural parameters implies a pair of low-frequency slopes. + +```{code-cell} ipython3 +h_draws = np.array([h_zero(unpack(v)) for v in kept]) + +fig, axes = plt.subplots(1, 2, figsize=(11, 3.8)) +for ax, j, ttl in zip(axes, [0, 1], + [r'$h_{\pi,\Delta m}(0)$', r'$h_{R,\Delta m}(0)$']): + ax.hist(h_draws[:, j], bins=60, color='C0', alpha=0.8, density=True) + ax.axvline(1, color='0.3', ls='--', lw=1.5) + ax.set_title(f'{ttl} mean {h_draws[:, j].mean():.3f}, ' + f'90% [{np.percentile(h_draws[:, j], 5):.3f}, ' + f'{np.percentile(h_draws[:, j], 95):.3f}]') +fig.suptitle('Posterior of the low-frequency slopes implied by the model') +plt.tight_layout() +plt.show() +``` + +The estimated model reproduces Lucas's two illustrations over a sample much like his. + +Compare these with the VAR estimates for 1960-1983 that we computed earlier, and with the paper's reported posterior means of 1.0068 and 0.8163. + +## How monetary policy moves the low-frequency slopes + +Now for the paper's main experiment. + +We lock every structural parameter at its posterior mean, vary only the two policy coefficients, and recompute $h(0)$ at each point. + +```{code-cell} ipython3 +p_bar = unpack(kept.mean(0)) + + +def h_grid(p, rule, k1, k2, g1_vals, g2_vals): + H1 = np.full((len(g2_vals), len(g1_vals)), np.nan) + H2 = np.full_like(H1, np.nan) + for i, b in enumerate(g2_vals): + for j, a in enumerate(g1_vals): + H1[i, j], H2[i, j] = h_zero(dict(p, **{k1: a, k2: b}), rule) + return H1, H2 + + +phi_pi_grid = np.linspace(-3, 1, 49) +phi_x_grid = np.linspace(-1, 0, 25) +Hpi, HR = h_grid(p_bar, 'money', 'phi_pi', 'phi_x', phi_pi_grid, phi_x_grid) +``` + +```{code-cell} ipython3 +def contour_panel(g1_vals, g2_vals, H1, H2, xlab, ylab, suptitle, scatter=None): + fig, axes = plt.subplots(1, 2, figsize=(12, 4.4), sharey=True) + for ax, Hm, ttl in zip(axes, [H1, H2], + [r'$h_{\pi,\Delta m}(0)$', r'$h_{R,\Delta m}(0)$']): + cs = ax.contourf(g1_vals, g2_vals, Hm, levels=14, cmap='viridis') + ax.contour(g1_vals, g2_vals, Hm, levels=[0.2, 1.0], + colors=['white', 'red'], linewidths=1.8) + ax.contourf(g1_vals, g2_vals, np.isnan(Hm).astype(float), + levels=[0.5, 1.5], colors=['0.85']) + if scatter is not None: + ax.scatter(*scatter, s=1.5, color='k', alpha=0.25) + fig.colorbar(cs, ax=ax) + ax.set(title=ttl, xlabel=xlab) + axes[0].set_ylabel(ylab) + fig.suptitle(suptitle) + plt.tight_layout() + plt.show() + + +contour_panel(phi_pi_grid, phi_x_grid, Hpi, HR, r'$\phi_\pi$', r'$\phi_x$', + 'Low-frequency slopes under the money growth rule\n' + '(white contour: 0.2, red contour: 1.0, dots: posterior draws)', + scatter=(kept[:, FREE.index('phi_pi')], + kept[:, FREE.index('phi_x')])) +``` + +Three things come through, and they are the paper's three findings. + +The low-frequency slopes are *not* policy invariant: moving $(\phi_\pi, \phi_x)$ alone sweeps $h_{\pi,\Delta m}(0)$ across most of the range that the rolling VARs displayed. + +The cloud of posterior draws sits in the region where both slopes are near one, which is where the 1960-1983 data put us. + +And a much more aggressive anti-inflation stance, meaning a large *negative* $\phi_\pi$ so that money growth is cut hard when inflation rises, moves both slopes far below one, with the smallest values in the corner where $\phi_x$ is also near zero. + +The red contour marks $h(0) = 1$ and, as {prf:ref}`ss_prop` requires, it passes exactly through $\phi_\pi = \phi_x = 0$. + +Any gray region would be one with no unique covariance stationary equilibrium, where $h(0)$ does not exist. + +Here is where the slopes bottom out, a slice through the grid at the posterior mean of $\phi_x$, and the fraction of the grid the stationarity test removes. + +```{code-cell} ipython3 +print(f'no stationary equilibrium at {np.isnan(Hpi).mean():.0%} of grid points') +for name, Hm in [('h_pi', Hpi), ('h_R', HR)]: + i, j = np.unravel_index(np.nanargmin(Hm), Hm.shape) + print(f' smallest {name:5s} = {Hm[i, j]:6.3f} at ' + f'phi_pi = {phi_pi_grid[j]:5.2f}, phi_x = {phi_x_grid[i]:5.2f}') + +row = np.argmin(np.abs(phi_x_grid - p_bar['phi_x'])) +print(f'\nslice at phi_x = {phi_x_grid[row]:.3f}') +for j in range(0, len(phi_pi_grid), 6): + print(f' phi_pi = {phi_pi_grid[j]:6.2f}: ' + f'h_pi = {Hpi[row, j]:7.3f}, h_R = {HR[row, j]:7.3f}') +``` + +### The same exercise with a Taylor rule + +The paper repeats the experiment with the interest rate rule {eq}`eq:ss_trule`. + +We keep the smoothing coefficient and the policy shock standard deviation at their money rule estimates and vary $(\psi_\pi, \psi_x)$. + +Where the Taylor rule leaves the equilibrium indeterminate we simply record it, shading those regions in gray; the paper instead selects the orthogonality solution of {cite:t}`LubikSchorfheide2004`. + +```{code-cell} ipython3 +p_taylor = dict(p_bar, rho_r=p_bar['rho_m'], psi_pi=1.5, psi_x=0.5) +psi_pi_grid = np.linspace(0, 3, 37) +psi_x_grid = np.linspace(0, 1, 21) +Tpi, TR = h_grid(p_taylor, 'taylor', 'psi_pi', 'psi_x', + psi_pi_grid, psi_x_grid) + +print(f'indeterminate or nonexistent at {np.isnan(Tpi).mean():.0%} of grid points') +print(f'h_pi ranges over [{np.nanmin(Tpi):.2f}, {np.nanmax(Tpi):.2f}]') +print(f'h_R ranges over [{np.nanmin(TR):.2f}, {np.nanmax(TR):.2f}]') + +contour_panel(psi_pi_grid, psi_x_grid, Tpi, TR, r'$\psi_\pi$', r'$\psi_x$', + 'Low-frequency slopes under a Taylor rule\n' + '(gray: no unique stable equilibrium)') +``` + +The dependence on policy survives, but it is much weaker. + +Over the determinate region the Taylor rule cannot drive either slope anywhere near zero. + +This is the paper's finding that the interest rate rule is "less able to replicate the estimates." + +The decomposition {eq}`eq:ss_decomp` supplies a reason. + +Driving $h_{\pi,\Delta m}(0)$ toward zero requires making the ratio $S_{z,\pi+z}(0)/S_{\pi+z}(0)$ close to one, which in practice means squeezing the zero-frequency power of inflation until $\pi_t + z_t$ is dominated by technology. + +A money growth rule with a large negative $\phi_\pi$ does exactly that, because the central bank is choosing the left side of {eq}`eq:ss_qtm` directly. + +A Taylor rule instead sets the interest rate and lets money growth be whatever money demand implies, so it never gets much purchase on low-frequency inflation relative to money growth, and $h_{\pi,\Delta m}(0)$ stays near one. + +## Estimating with Hamiltonian Monte Carlo + +The random walk sampler that produced our posterior is the traditional workhorse of +DSGE estimation. + +It is also wasteful, and we can measure how wasteful. + +This section takes the Sargent-Surico model as a guinea pig for a modern +alternative. + +### What the random walk delivers + +A correlated chain of length $N$ is worth fewer than $N$ independent draws, and the +**effective sample size** says how many. + +```{code-cell} ipython3 +import arviz as az + +rwmh_idata = az.from_dict({n: kept[:, j][None, :] for j, n in enumerate(FREE)}) +rwmh_summary = az.summary(rwmh_idata, var_names=FREE) +print(rwmh_summary[['mean', 'sd', 'ess_bulk']].to_string()) +print(f'\n{len(kept)} retained draws out of {N_DRAWS}, in {rwmh_seconds:.0f} seconds') +print(f'smallest effective sample size = {rwmh_summary["ess_bulk"].min():.0f}') +``` + +Tens of thousands of iterations buy us on the order of a hundred independent draws +for the worst-mixing parameter. + +The reason is visible in the proposal. + +A random walk perturbs all eighteen parameters in a direction that knows nothing +about the local shape of the posterior, so the step size has to respect the most +tightly constrained direction while the loosest direction is explored at that same +crawl. + +Our posterior is badly scaled. + +```{code-cell} ipython3 +w = np.linalg.eigvalsh(H) +print(f'condition number of the posterior Hessian {w.max() / w.min():10.3g}') +print(f'ratio of longest to shortest posterior scale {np.sqrt(w.max() / w.min()):10.0f}') +``` + +### Why Hamiltonian methods are awkward for DSGE models + +**Hamiltonian Monte Carlo** {cite}`DuaneEtAl1987,Neal2011` treats the parameter +vector as the position of a particle, gives it a random momentum, and simulates +Hamiltonian dynamics in which the negative log posterior plays the role of potential +energy. + +Trajectories follow the contours of the distribution instead of blundering across +them, so a proposal can travel a long way and still be accepted. + +The **No-U-Turn sampler** {cite}`HoffmanGelman2014` removes the need to tune the +trajectory length by extending each trajectory until it starts to double back on +itself. + +{cite:t}`Betancourt2017` gives a conceptual introduction. + +The price of admission is the gradient $\nabla_\theta \log p(\theta \mid Y)$. + +That is where DSGE models become awkward. + +Our likelihood runs through `gensys`, whose first act is an *ordered* generalized +Schur decomposition: it computes generalized eigenvalues, **sorts** them by modulus, +and reorders the decomposition to match. + +Sorting is a discrete operation, and no automatic differentiation library propagates +derivatives through it. + +This, rather than any statistical objection, is the reason Hamiltonian methods +remain uncommon in DSGE estimation. + +### A differentiable solver + +The obstacle is the *algorithm*, not the model, so we change algorithms. + +Write the equilibrium conditions as + +$$ +F_A\, \mathbb{E}_t y_{t+1} + F_B\, y_t + F_C\, y_{t-1} + F_E\, \varepsilon_t = 0 , +$$ (eq:ss_structural) + +where + +$$ +y_t = [\pi_t,\ x_t,\ \Delta m_t,\ R_t,\ e_t,\ a_t,\ \chi_t,\ z_t,\ u_t]' +$$ + +and $u_t = \varepsilon_{mt}$ carries the money rule shock. + +The auxiliary expectation variables that `gensys` needs are not required here. + +Conjecture a solution + +$$ +y_t = G\, y_{t-1} + \Theta\, \varepsilon_t . +$$ (eq:ss_conjecture) + +Then $\mathbb{E}_t y_{t+1} = G y_t$, and {eq}`eq:ss_structural` becomes + +$$ +(F_A G + F_B)\, y_t + F_C\, y_{t-1} + F_E\, \varepsilon_t = 0 . +$$ + +Solving for $y_t$ and matching coefficients with {eq}`eq:ss_conjecture` gives a fixed +point in $G$ alone, + +$$ +G = -(F_A G + F_B)^{-1} F_C , +\qquad +\Theta = -(F_A G + F_B)^{-1} F_E . +$$ (eq:ss_fixedpoint) + +Iterating {eq}`eq:ss_fixedpoint` from $G = 0$ involves nothing but matrix products +and linear solves, every one of them differentiable. + +```{code-cell} ipython3 +import jax +import jax.numpy as jnp +import numpyro +import numpyro.distributions as dist +from jax import lax +from numpyro.infer import MCMC, NUTS + +jax.config.update('jax_enable_x64', True) + +U, NJ = 8, 9 # y = [pi, x, dm, R, e, a, chi, z, u] + + +def structural(p): + """Return F_A, F_B, F_C, F_E of equation (SS-structural).""" + th, api, kap = p['theta'], p['alpha_pi'], p['kappa'] + ax, sig, xi, gam = p['alpha_x'], p['sigma'], p['xi'], p['gamma'] + rm, ra = p['rho_m'], p['rho_a'] + FA = jnp.zeros((NJ, NJ)); FB = jnp.zeros((NJ, NJ)) + FC = jnp.zeros((NJ, NJ)); FE = jnp.zeros((NJ, 5)) + + FA = FA.at[0, PI].set(-th * (1 - api)) # Phillips curve + FB = FB.at[0, PI].set(1).at[0, X].set(-kap).at[0, E].set(1 / TAU) + FC = FC.at[0, PI].set(-th * api) + + FA = FA.at[1, X].set(-(1 - ax)).at[1, PI].set(-sig) # IS curve + FB = (FB.at[1, X].set(1).at[1, R].set(sig) + .at[1, A_].set(-sig * (1 - xi) * (1 - ra))) + FC = FC.at[1, X].set(-ax) + + FB = FB.at[2, DM].set(1).at[2, PI].set(-1).at[2, Z].set(-1) # money demand + FB = (FB.at[2, X].set(-1 / (sig * gam)).at[2, R].set(1 / gam) + .at[2, CH].set(-1 / gam).at[2, A_].set(1 / gam)) + FC = (FC.at[2, X].set(1 / (sig * gam)).at[2, R].set(-1 / gam) + .at[2, CH].set(1 / gam).at[2, A_].set(-1 / gam)) + + FB = (FB.at[3, DM].set(1).at[3, U].set(-1) # policy rule + .at[3, PI].set(-(1 - rm) * p['phi_pi']) + .at[3, X].set(-(1 - rm) * p['phi_x'])) + FC = FC.at[3, DM].set(-rm) + + for row, (v, rho, k) in enumerate([(E, p['rho_e'], 0), (A_, ra, 1), + (CH, p['rho_chi'], 2), (Z, 0.0, 3)], 4): + FB = FB.at[row, v].set(1) + FC = FC.at[row, v].set(-rho) + FE = FE.at[row, k].set(-1) + FB = FB.at[8, U].set(1) + FE = FE.at[8, 4].set(-1) + return FA, FB, FC, FE + + +def solve_fixed_point(p, n_iter=60): + """Iterate (SS-fixedpoint) to convergence; differentiable throughout.""" + FA, FB, FC, FE = structural(p) + + def step(G, _): + return -jnp.linalg.solve(FA @ G + FB, FC), None + + G, _ = lax.scan(step, jnp.zeros((NJ, NJ)), None, length=n_iter) + return G, -jnp.linalg.solve(FA @ G + FB, FE) +``` + +Before trusting it we check it against `gensys` at the paper's posterior means. + +The two algorithms share no code, so agreement to machine precision is a real test. + +```{code-cell} ipython3 +p_check = {n: float(PAPER[n]) for n in FREE} +G_fp, Theta_fp = solve_fixed_point(p_check) + +A_gen, B_gen, C_gen, _ = state_space(dict(p_check, tau=TAU)) +idx = [PI, X, DM, R, E, A_, CH, Z] +sd_check = np.array([p_check[n] for n in + ['sig_e', 'sig_a', 'sig_chi', 'sig_z', 'sig_m']]) + +print('fixed point versus gensys') +print(f' transition matrix max abs difference ' + f'{np.abs(np.asarray(G_fp)[np.ix_(idx, idx)] - A_gen[:10, :10][np.ix_(idx, idx)]).max():.2e}') +print(f' shock loadings max abs difference ' + f'{np.abs(np.asarray(Theta_fp * sd_check)[idx] - B_gen[:10][idx]).max():.2e}') +``` + +```{note} +The fixed point converges to the stable solution when there is one, but unlike +`gensys` it returns no existence and uniqueness flags. + +We therefore keep `gensys` for the determinacy map of the previous section and use +the fixed point only where we need derivatives. +``` + +### The likelihood in JAX + +The Kalman filter carries over unchanged, written with `lax.scan` so that it too +can be differentiated. + +The only substantive change is the initial covariance: instead of calling a +Lyapunov solver we use $\operatorname{vec}(P) = (I - A \otimes A)^{-1}\operatorname{vec}(Q)$, +which is one linear solve. + +```{code-cell} ipython3 +def state_space_jax(p): + G, Theta = solve_fixed_point(p) + sd = jnp.array([p['sig_e'], p['sig_a'], p['sig_chi'], p['sig_z'], p['sig_m']]) + n = NJ + 2 # append x(-1) and a(-1) + A = jnp.zeros((n, n)).at[:NJ, :NJ].set(G) + A = A.at[NJ, X].set(1).at[NJ + 1, A_].set(1) + B = jnp.zeros((n, 5)).at[:NJ].set(Theta * sd) + C = jnp.zeros((4, n)) + C = C.at[0, DM].set(1).at[1, PI].set(1).at[2, R].set(1) + C = C.at[3, X].set(1).at[3, NJ].set(-1) + C = C.at[3, A_].set(p['xi']).at[3, NJ + 1].set(-p['xi']) + C = C.at[3, Z].set(1) + return A, B, C + + +def loglik_jax(p, Y): + A, B, C = state_space_jax(p) + n = A.shape[0] + Q = B @ B.T + P0 = jnp.linalg.solve(jnp.eye(n * n) - jnp.kron(A, A), + Q.reshape(-1)).reshape(n, n) + const = Y.shape[1] * jnp.log(2 * jnp.pi) + + def step(carry, y): + s, P = carry + s, P = A @ s, A @ P @ A.T + Q + v = y - C @ s + PCt = P @ C.T + F = C @ PCt + L = jnp.linalg.cholesky(F) + u = jax.scipy.linalg.solve_triangular(L, v, lower=True) + ll = -0.5 * (const + 2 * jnp.sum(jnp.log(jnp.diag(L))) + u @ u) + K = jnp.linalg.solve(F, PCt.T).T + return (s + K @ v, P - K @ PCt.T), ll + + _, lls = lax.scan(step, (jnp.zeros(n), P0), Y) + return jnp.sum(lls) +``` + +Two checks: the new likelihood must agree with the one we used for the random walk, +and its gradient must agree with finite differences. + +```{code-cell} ipython3 +Y_jax = jnp.asarray(Y_est) +print(f'log likelihood, JAX {float(loglik_jax(p_check, Y_jax)):.8f}') +print(f'log likelihood, NumPy {loglik(dict(p_check, tau=TAU), Y_est):.8f}') + +grad_ll = jax.jit(jax.grad(lambda q: loglik_jax(q, Y_jax))) +g = grad_ll(p_check) + +print('\n autodiff finite difference') +for n in ['phi_pi', 'kappa', 'rho_m', 'sig_z']: + step_n = 1e-5 + up = dict(p_check); up[n] = p_check[n] + step_n + dn = dict(p_check); dn[n] = p_check[n] - step_n + fd = (loglik_jax(up, Y_jax) - loglik_jax(dn, Y_jax)) / (2 * step_n) + print(f' {n:8s} {float(g[n]):12.5f} {float(fd):17.5f}') +``` + +```{code-cell} ipython3 +t0 = time.time() +for _ in range(20): + gg = grad_ll(p_check) +jax.block_until_ready(gg) +grad_ms = 1000 * (time.time() - t0) / 20 + +t0 = time.time() +for _ in range(20): + loglik(dict(p_check, tau=TAU), Y_est) +loglik_ms = 1000 * (time.time() - t0) / 20 + +print(f'one NumPy log likelihood {loglik_ms:6.2f} ms') +print(f'one JAX gradient (18 partials) {grad_ms:6.2f} ms') +``` + +A full gradient costs about what a single likelihood evaluation costs, which is the +whole point of reverse-mode differentiation. + +A finite-difference gradient would need at least nineteen likelihood evaluations. + +### Sampling with NUTS + +We hand the same priors to NumPyro {cite}`PhanEtAl2019`, which supplies NUTS and +handles the transformations to unconstrained space that Hamiltonian dynamics +require. + +```{code-cell} ipython3 +def beta_np(m, s): + nu = m * (1 - m) / s ** 2 - 1 + return dist.Beta(m * nu, (1 - m) * nu) + + +def gamma_np(m, s): + return dist.Gamma(m ** 2 / s ** 2, m / s ** 2) + + +def invgamma_np(m, s): + a = m ** 2 / s ** 2 + 2 + return dist.InverseGamma(a, m * (a - 1)) + + +PRIORS_NP = { + 'theta': beta_np(0.99, 0.005), 'alpha_pi': beta_np(0.5, 0.2), + 'alpha_x': beta_np(0.5, 0.2), 'rho_m': beta_np(0.5, 0.05), + 'rho_e': beta_np(0.5, 0.1), 'rho_a': beta_np(0.5, 0.1), + 'rho_chi': beta_np(0.5, 0.1), + 'kappa': gamma_np(0.3, 0.1), 'sigma': gamma_np(0.1, 0.05), + 'xi': gamma_np(2.0, 1.0), 'gamma': gamma_np(4.0, 1.0), + 'phi_pi': dist.Normal(0, 0.5), 'phi_x': dist.Normal(0, 0.5), + 'sig_e': invgamma_np(0.3, 1.0), 'sig_a': invgamma_np(0.3, 1.0), + 'sig_chi': invgamma_np(0.3, 1.0), 'sig_z': invgamma_np(0.3, 1.0), + 'sig_m': invgamma_np(0.3, 1.0)} + + +def ss_model(Y): + p = {n: numpyro.sample(n, PRIORS_NP[n]) for n in FREE} + numpyro.factor('loglik', loglik_jax(p, Y)) +``` + +Three settings matter. + +We ask for a **dense mass matrix**, because the condition number reported above says +the posterior has correlations that a diagonal preconditioner cannot absorb. + +We cap the trajectory length, since without a cap NUTS spends most of its time on very +long trajectories in the flattest directions. + +And we run **four chains** rather than one, all started from the posterior mode. + +That last choice is the one that earns its keep. + +```{code-cell} ipython3 +kernel = NUTS(ss_model, target_accept_prob=0.8, dense_mass=True, + max_tree_depth=8, + init_strategy=numpyro.infer.init_to_value( + values={n: float(v) for n, v in zip(FREE, v_mode)})) +mcmc = MCMC(kernel, num_warmup=400, num_samples=400, num_chains=4, + chain_method='sequential', progress_bar=False) + +t0 = time.time() +mcmc.run(jax.random.PRNGKey(1), Y_jax, extra_fields=('num_steps', 'diverging')) +jax.block_until_ready(mcmc.get_samples()) +nuts_seconds = time.time() - t0 + +extra = mcmc.get_extra_fields() +print(f'{nuts_seconds:.0f} seconds for 4 chains of 400 draws') +print(f'mean leapfrog steps per iteration ' + f'{np.asarray(extra["num_steps"]).mean():.0f}') +print(f'divergences ' + f'{int(np.asarray(extra["diverging"]).sum())}') +``` + +```{code-cell} ipython3 +nuts_idata = az.from_numpyro(mcmc) +nuts_summary = az.summary(nuts_idata, var_names=FREE) +nuts_kept = np.column_stack([np.asarray(mcmc.get_samples()[n]) for n in FREE]) +nuts_summary[['mean', 'sd', 'hdi_3%', 'hdi_97%', 'ess_bulk', 'r_hat']] +``` + +### A warning sign in the diagnostics + +Most parameters look excellent, with $\hat R$ at one and effective sample sizes in the +hundreds or thousands. + +A few do not, and it is worth asking which. + +```{code-cell} ipython3 +worst = nuts_summary['ess_bulk'].nsmallest(3).index.tolist() +print('weakest mixing:') +print(nuts_summary.loc[worst, ['mean', 'sd', 'ess_bulk', 'r_hat']]) + +by_chain = mcmc.get_samples(group_by_chain=True) +print('\nper-chain posterior means') +print(f'{"":10s}' + ''.join(f'{"chain " + str(c):>10s}' for c in range(4))) +for n in worst: + v = np.asarray(by_chain[n]) + print(f'{n:10s}' + ''.join(f'{v[c].mean():10.4f}' for c in range(4))) +``` + +The chains started from a common point, so any disagreement between them means some +of them have wandered somewhere the others have not. + +The identity of the badly behaved parameters is a clue. + +Inflation in the data is persistent, and equation {eq}`eq:ss_nkpc` offers two ways to +deliver that persistence. + +**Intrinsic persistence** comes from indexation to past inflation, a large +$\alpha_\pi$, with the markup shock left transitory. + +**Inherited persistence** comes from a persistent markup shock, a large $\rho_e$, with +little indexation. + +The parameters that mix worst are exactly the ones that distinguish these two stories. + +That is a reason to suspect the posterior has more than one mode, and it tells us +where to look for the second one. + +### Two modes + +Sampler diagnostics are a noisy instrument, so we settle the question with the +optimizer instead. + +Our mode search in the estimation section started from the paper's posterior means and +climbed to a high-indexation mode. + +We now start the same optimizer from the opposite corner. + +```{code-cell} ipython3 +v_start_a = v_mode.copy() +for n, val in [('alpha_pi', 0.12), ('rho_e', 0.84), ('sig_e', 0.35)]: + v_start_a[FREE.index(n)] = val + +res_a = optimize.minimize(neg_log_post, v_start_a, method='Powell', + options=dict(maxiter=20000, maxfev=20000)) +res_a = optimize.minimize(neg_log_post, res_a.x, method='L-BFGS-B', + bounds=[SUPPORT[n] for n in FREE]) +v_mode_a = res_a.x + +print(f'log posterior, low-indexation mode {-res_a.fun:10.3f}') +print(f'log posterior, high-indexation mode {log_post(v_mode, Y_est):10.3f}') +print('\n low-index high-index paper') +for n in ['alpha_pi', 'rho_e', 'sig_e', 'kappa', 'phi_pi', 'phi_x']: + i = FREE.index(n) + print(f' {n:9s}{v_mode_a[i]:10.4f}{v_mode[i]:13.4f}{PAPER[n]:11.4f}') +``` + +There are two distinct modes, and the low-indexation one has the **higher** +posterior density. + +The mode we found earlier, the one nearest the paper's published estimates and the one +our random walk explored for thirty thousand draws, is a *local* mode. + +The pooled NUTS draws show the same thing from the sampling side. + +```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Marginal posteriors for the parameters that separate the two modes + name: fig-ss-bimodal +--- +fig, axes = plt.subplots(1, 3, figsize=(13, 3.8)) +for ax, n, ttl in zip(axes, ['alpha_pi', 'rho_e', 'sig_e'], + [r'$\alpha_\pi$: indexation', + r'$\rho_e$: markup persistence', + r'$\sigma_e$: markup volatility']): + i = FREE.index(n) + ax.hist(nuts_kept[:, i], bins=60, density=True, color='C0', alpha=0.8) + ax.axvline(v_mode_a[i], color='C2', lw=2, label='low-index mode') + ax.axvline(v_mode[i], color='C3', lw=2, label='high-index mode') + ax.axvline(PAPER[n], color='k', ls='--', lw=1.2, label='paper') + ax.set(title=ttl, xlabel=n) +axes[0].legend(fontsize=8) +fig.suptitle('The marginals spread across both modes rather than concentrating on one') +fig.tight_layout() +plt.show() +``` + +The marginals are wide and lumpy rather than the tidy bells that every other parameter +produces. + +The chains do move between the two regions, but slowly, and that slow movement is +precisely what the poor effective sample sizes were reporting. + +Thirty thousand random walk draws never made the trip at all. + +The economics behind the two modes is a familiar identification problem. + +The data over 1960-1983 cannot tell us whether inflation is persistent because prices +are indexed to past inflation or because the shocks pushing inflation around are +themselves persistent. + +Both stories fit, and the published estimates describe the one the data like slightly +less. + +### Does the bimodality change the economics? + +This is the question that matters for the rest of the lecture. + +```{code-cell} ipython3 +h_mode_a = h_zero(unpack(v_mode_a)) +h_mode_b = h_zero(unpack(v_mode)) +print(f'h(0) at the low-indexation mode {h_mode_a[0]:.4f}, {h_mode_a[1]:.4f}') +print(f'h(0) at the high-indexation mode {h_mode_b[0]:.4f}, {h_mode_b[1]:.4f}') + +h_nuts = np.array([h_zero(unpack(v)) for v in nuts_kept]) +low = nuts_kept[:, FREE.index('alpha_pi')] < 0.5 +print(f'\n{low.sum()} of {len(low)} NUTS draws sit on the low-indexation side') +for tag, m in [('low-indexation draws ', low), ('high-indexation draws', ~low)]: + if m.sum() < 20: + continue + print(f'{tag} h_pi {h_nuts[m, 0].mean():.4f} ' + f'[{np.percentile(h_nuts[m, 0], 5):.4f}, {np.percentile(h_nuts[m, 0], 95):.4f}]' + f' h_R {h_nuts[m, 1].mean():.4f} ' + f'[{np.percentile(h_nuts[m, 1], 5):.4f}, {np.percentile(h_nuts[m, 1], 95):.4f}]') +``` + +```{code-cell} ipython3 +fig, axes = plt.subplots(1, 2, figsize=(11, 3.8)) +for ax, j, ttl in zip(axes, [0, 1], + [r'$h_{\pi,\Delta m}(0)$', r'$h_{R,\Delta m}(0)$']): + ax.hist(h_draws[:, j], bins=50, density=True, alpha=0.55, + color='C0', label='RWMH') + ax.hist(h_nuts[:, j], bins=50, density=True, alpha=0.55, + color='C3', label='NUTS') + ax.axvline(1, color='0.3', ls='--', lw=1.5) + ax.set_title(ttl) + ax.legend(fontsize=9) +fig.suptitle('The low-frequency slopes are robust across samplers and modes') +plt.tight_layout() +plt.show() +``` + +The reassuring answer is that it does not. + +The two modes disagree sharply about *why* inflation is persistent and far less about +the low-frequency slopes that this lecture is about. + +Whether inflation persistence is intrinsic or inherited, $h_{\pi,\Delta m}(0)$ stays +between about four fifths and one, and $h_{R,\Delta m}(0)$ sits near four fifths in +both cases. + +The paper's substantive conclusions survive the discovery that its parameter estimates +describe a local mode. + +### Comparing the samplers + +With that settled we can return to efficiency. + +Effective sample sizes for the three parameters that separate the modes are not +comparable across the two samplers, because for NUTS they partly measure movement +*between* modes while for the random walk they measure movement within one. + +We therefore compare medians rather than minima. + +```{code-cell} ipython3 +compare = pd.DataFrame({ + 'RWMH mean': rwmh_summary['mean'], 'NUTS mean': nuts_summary['mean'], + 'RWMH ESS': rwmh_summary['ess_bulk'], 'NUTS ESS': nuts_summary['ess_bulk'], + 'NUTS r_hat': nuts_summary['r_hat']}) +compare['ESS ratio'] = compare['NUTS ESS'] / compare['RWMH ESS'] +compare.round(3) +``` + +```{code-cell} ipython3 +unimodal = ~nuts_summary.index.isin(worst) +r_med = rwmh_summary['ess_bulk'][unimodal].median() +n_med = nuts_summary['ess_bulk'][unimodal].median() + +print(f'{"":26s}{"RWMH":>12s}{"NUTS":>12s}') +print(f'{"draws kept":26s}{len(kept):12d}{len(nuts_kept):12d}') +print(f'{"seconds":26s}{rwmh_seconds:12.0f}{nuts_seconds:12.0f}') +print(f'{"median ESS":26s}{r_med:12.0f}{n_med:12.0f}') +print(f'{"median ESS per 1000 draws":26s}{1000 * r_med / len(kept):12.1f}' + f'{1000 * n_med / len(nuts_kept):12.1f}') +print(f'{"median ESS per second":26s}{r_med / rwmh_seconds:12.2f}' + f'{n_med / nuts_seconds:12.2f}') +print(f'\n{"speed-up in ESS per second":30s}' + f'{(n_med / nuts_seconds) / (r_med / rwmh_seconds):6.1f}x') +``` + +Per draw the gap is very large, because each NUTS draw ends a trajectory that has +travelled across the posterior rather than taking one small random step. + +Cost claws some of that back, since a NUTS iteration needs many gradient evaluations +where a random walk iteration needs one likelihood evaluation. + +Effective draws per second is the ratio that matters, and it still favours NUTS by a +wide margin. + +```{code-cell} ipython3 +:tags: [hide-input] + +fig, axes = plt.subplots(1, 2, figsize=(12, 4.5)) +order = np.argsort(rwmh_summary['ess_bulk'].values) +pos = np.arange(len(FREE)) +axes[0].barh(pos - 0.2, rwmh_summary['ess_bulk'].values[order], 0.4, + label='RWMH', color='C0') +axes[0].barh(pos + 0.2, nuts_summary['ess_bulk'].values[order], 0.4, + label='NUTS', color='C3') +axes[0].set(yticks=pos, yticklabels=[FREE[i] for i in order], + xlabel='effective sample size', title='ESS by parameter') +axes[0].legend() + +j = FREE.index('rho_a') +axes[1].plot(np.linspace(0, 1, len(kept)), kept[:, j], lw=0.5, color='C0', + alpha=0.8, label=f'RWMH ({len(kept)} draws)') +axes[1].plot(np.linspace(0, 1, len(nuts_kept)), nuts_kept[:, j], lw=0.5, + color='C3', alpha=0.8, label=f'NUTS ({len(nuts_kept)} draws)') +axes[1].set(xlabel='fraction of the run', ylabel=r'$\rho_a$', + title=r'traces for $\rho_a$, a parameter both samplers agree on') +axes[1].legend() +fig.tight_layout() +plt.show() +``` + +### What to take away + +The obstacle to Hamiltonian methods in DSGE estimation is the model solver, not the +statistics. + +Replacing an eigenvalue-sorting solver by a fixed point that uses only linear algebra +buys exact gradients for about the cost of one extra likelihood evaluation, and that +is enough to put NUTS within reach. + +The payoff was not only speed. + +Cheap chains made it cheap to run several of them and inspect their diagnostics, and +the parameters that mixed worst pointed straight at a second and better mode that +thirty thousand random walk draws had never visited. + +Two caveats are worth stating. + +First, the fixed point returns something meaningless where no stable solution exists, +so a sampler that wanders outside the determinacy region will produce nonsense rather +than a rejection; the posterior here sits well inside that region, but a flatter one +might not. + +Second, the efficiency gain reported above is specific to this posterior, and a +well-scaled posterior in fewer dimensions would narrow the gap considerably. + +## Concluding remarks + +{cite:t}`Lucas1980` was careful to say that his two unit slopes should be expected to hold under some monetary policies and to break down under others. + +Extending his sample by half a century confirms the breakdowns. + +Estimating a small structural model on a sample like his, and then moving only the monetary policy rule, produces slopes that range over most of what the data display. + +The lesson is the one Lucas, {cite:t}`Sargent1971`, and {cite:t}`Whiteman1984` all drew, namely that low-frequency regression coefficients are not structural. + +Working through the model from scratch adds something to the paper's account. + +Lucas's first illustration, the unit slope of inflation on money growth, is an *identity* in this model whenever money growth is econometrically exogenous, because money demand is a quantity equation whose remaining terms are all first differences. + +His second illustration, the unit slope of the interest rate on money growth, has no such support and moves around freely. + +So the two illustrations, which look symmetric in a scatter plot, are quite different objects from the perspective of a structural model. + +That asymmetry also marks the limit of the exercise, as {ref}`ss_ex2` brings out. + +Because $h_{\pi,\Delta m}(0)$ is anchored at one by money demand, driving it to the near-zero values that the post-1984 data display takes a money growth rule far more anti-inflationary than the one those same data select. + +On the computational side, the model turned out to be a useful test bed for Hamiltonian Monte Carlo. + +The barrier to using it on DSGE models is not statistical but algorithmic: the standard solvers sort eigenvalues, and sorting has no derivative. + +Swapping in a fixed-point solver that uses only linear algebra restores exact gradients, and the sampler that becomes available is far more efficient per unit of computing time on a posterior as badly scaled as this one. + +The larger dividend was a diagnostic one. + +Cheap chains and their convergence diagnostics led us to a second mode, and the posterior turned out to be bimodal: the data cannot choose between intrinsic and inherited inflation persistence, and the mode nearest the published estimates is the *lower* of the two. + +That the two modes nonetheless agree about $h_{\pi,\Delta m}(0)$ and $h_{R,\Delta m}(0)$ is a piece of good news for the paper that only the more thorough sampler could deliver. + +## Exercises + +```{exercise} +:label: ss_ex1 + +The paper takes technology growth to be serially uncorrelated. + +Replace $z_t = \varepsilon_{zt}$ by $z_t = \rho_z z_{t-1} + \varepsilon_{zt}$ and ask two questions. + +Does {prf:ref}`ss_prop` survive, and does persistent technology growth offer a route to breaking the unit slope that has nothing to do with monetary policy? +``` + +```{solution-start} ss_ex1 +:class: dropdown +``` + +The technology block is row 7 of the canonical form, where `g1[7, Z]` is currently zero. + +```{code-cell} ipython3 +def canonical_rho_z(p, rule='money'): + g0, g1, psi, pie = canonical(p, rule) + g1[7, Z] = p['rho_z'] + return g0, g1, psi, pie + + +def h_zero_rho_z(p): + G1, impact, eu = gensys(*canonical_rho_z(p)) + if eu != (1, 1): + return np.nan, np.nan + sd = np.diag([p['sig_e'], p['sig_a'], p['sig_chi'], p['sig_z'], p['sig_m']]) + n = NY + 2 + A, B = np.zeros((n, n)), np.zeros((n, 5)) + A[:NY, :NY], A[NY, X], A[NY + 1, A_] = G1, 1, 1 + B[:NY] = impact @ sd + C = np.zeros((4, n)) + C[0, DM], C[1, PI], C[2, R], C[3, Z] = 1, 1, 1, 1 + C[3, X], C[3, NY] = 1, -1 + C[3, A_], C[3, NY + 1] = p['xi'], -p['xi'] + if np.max(np.abs(np.linalg.eigvals(A))) > 1 - 1e-9: + return np.nan, np.nan + return (h_zero_from_state_space(A, B, C, 1, 0), + h_zero_from_state_space(A, B, C, 2, 0)) + + +print('exogenous money growth, phi_pi = phi_x = 0') +for rho_z in [0.0, 0.3, 0.6, 0.9]: + q = dict(p_bar, phi_pi=0.0, phi_x=0.0, rho_z=rho_z) + print(f' rho_z = {rho_z:.1f}: h_pi = {h_zero_rho_z(q)[0]:.10f}, ' + f'h_R = {h_zero_rho_z(q)[1]:.4f}') + +print('\nwith policy feedback, phi_pi = -2, phi_x = -0.5') +for rho_z in [0.0, 0.3, 0.6, 0.9]: + q = dict(p_bar, phi_pi=-2.0, phi_x=-0.5, rho_z=rho_z) + print(f' rho_z = {rho_z:.1f}: h_pi = {h_zero_rho_z(q)[0]:.4f}, ' + f'h_R = {h_zero_rho_z(q)[1]:.4f}') +``` + +{prf:ref}`ss_prop` survives untouched, because its proof used only that technology is orthogonal to the policy shock, never how persistent it is. + +Under exogenous money growth *both* slopes are invariant to $\rho_z$, and for a reason worth noticing: $z_t$ enters the model only contemporaneously, through money demand, so $\rho_z$ changes the dynamics of technology without changing anyone's response to a money shock. + +Once policy feeds back the picture is different, and persistent technology growth pulls $h_{\pi,\Delta m}(0)$ well away from one. + +So the answer to the second question is yes, and it points at a limitation of the paper's experiment: the technology process is held fixed while policy varies, but the two interact. + +```{solution-end} +``` + +```{exercise} +:label: ss_ex2 + +The lecture estimates the model on 1960:I-1983:IV. + +Find the posterior mode on 1984:I-2007:IV and compare the policy coefficients and the implied $h(0)$ with what we found for the earlier sample. + +The VAR put $h_{\pi,\Delta m}(0)$ near zero over 1984-2007, so does the structural model, re-estimated, agree? +``` + +```{solution-start} ss_ex2 +:class: dropdown +``` + +```{code-cell} ipython3 +mask2 = (data.index.year >= 1984) & (data.index.year <= 2007) +Y2 = data.loc[mask2, VARS].values +Y2 = Y2 - Y2.mean(0) + +res2 = optimize.minimize(lambda v: -log_post(v, Y2), v_mode, method='Powell', + options=dict(maxiter=20000, maxfev=20000)) +res2 = optimize.minimize(lambda v: -log_post(v, Y2), res2.x, method='L-BFGS-B', + bounds=[SUPPORT[n] for n in FREE]) +v2 = res2.x +print(' 1960-1983 1984-2007') +for n in ['phi_pi', 'phi_x', 'rho_m']: + j = FREE.index(n) + print(f' {n:9s} {v_mode[j]:11.4f} {v2[j]:12.4f}') +print(' %-9s %5.3f, %.3f %5.3f, %.3f' + % (('model h(0)',) + h_zero(unpack(v_mode)) + h_zero(unpack(v2)))) +for name, lo, hi in [('1960-1983', 1960, 1983), ('1984-2007', 1984, 2007)]: + m = (data.index.year >= lo) & (data.index.year <= hi) + med = np.nanmedian(bvar_h0(data.loc[m, VARS].values), axis=0) + print(f' VAR h(0) on {name}: {med[0]:.3f}, {med[1]:.3f}') +``` + +The re-estimated model does *not* agree with the VAR. + +The policy coefficients move, but the implied $h_{\pi,\Delta m}(0)$ stays close to one on both samples, while the VAR puts it near zero after 1984. + +{prf:ref}`ss_prop` explains why this had to be hard. + +Within this model $h_{\pi,\Delta m}(0)$ is anchored at one and is dragged away from it only by policy feedback, so reaching values near zero needs a large negative $\phi_\pi$, and the post-1984 data do not put the money growth rule anywhere near there. + +That is a genuine tension between the paper's Section II and its Section III: the contour plot shows a policy rule that *could* generate the low slopes, but the rule the later data actually select is not that rule. + +A full comparison would run the sampler on both samples rather than compare modes. + +```{solution-end} +``` + +```{exercise} +:label: ss_ex3 + +Figure 8 of the paper asks whether a fall in the variance of supply shocks, rather than a change in policy, could account for the decline in $h(0)$. + +Redo the money rule contour plot with $\sigma_e$ cut to one quarter of its posterior mean and report how much the picture moves. +``` + +```{solution-start} ss_ex3 +:class: dropdown +``` + +```{code-cell} ipython3 +p_small_e = dict(p_bar, sig_e=p_bar['sig_e'] / 4) +Hpi2, HR2 = h_grid(p_small_e, 'money', 'phi_pi', 'phi_x', + phi_pi_grid, phi_x_grid) + +print('baseline h_pi in [%.2f, %.2f], h_R in [%.2f, %.2f]' + % (np.nanmin(Hpi), np.nanmax(Hpi), np.nanmin(HR), np.nanmax(HR))) +print('small sig_e h_pi in [%.2f, %.2f], h_R in [%.2f, %.2f]' + % (np.nanmin(Hpi2), np.nanmax(Hpi2), np.nanmin(HR2), np.nanmax(HR2))) +print('largest change in h_pi across the grid: %.3f' + % np.nanmax(np.abs(Hpi2 - Hpi))) + +contour_panel(phi_pi_grid, phi_x_grid, Hpi2, HR2, r'$\phi_\pi$', r'$\phi_x$', + r'Low-frequency slopes with $\sigma_e$ cut to one quarter') +``` + +The range that $h_{\pi,\Delta m}(0)$ can reach is unchanged, and the mapping from policy to $h(0)$ keeps the same shape. + +The contours do move within that range, by up to a few tenths at some grid points, so the claim is not that shock variances are irrelevant. + +It is the weaker claim the paper makes, that a change in shock variances of this magnitude is not enough to substitute for a change in policy. + +```{solution-end} +``` diff --git a/lectures/var_subsets.md b/lectures/var_subsets.md new file mode 100644 index 000000000..f15ddcb41 --- /dev/null +++ b/lectures/var_subsets.md @@ -0,0 +1,1234 @@ +--- +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.16.7 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +(var_subsets)= +```{raw} jupyter +
+ + QuantEcon + +
+``` + +# Vector Autoregressions for Subsets of Variables + +```{index} single: Vector Autoregression; subsystems +``` + +```{index} single: Kalman Filter; and vector autoregressions +``` + +```{contents} Contents +:depth: 2 +``` + +In addition to what's in Anaconda, this lecture will need the following libraries: + +```{code-cell} ipython3 +:tags: [hide-output] + +!pip install quantecon +``` + +## Overview + +An economic model delivers a vector autoregression for a list of variables $Y_t$. + +An econometrician often observes only *some* of them. + +This lecture answers three questions about that situation. + +Given an $m$th order VAR for an $n \times 1$ vector $Y_t$ and a selector matrix +$S_y$ that extracts an $n_y \times 1$ subvector $y_t = S_y Y_t$: + +1. What vector autoregression does $y_t$ obey? +2. What is its moving average representation? +3. How are the innovations in the small system related to the innovations in the + large one? + +The answers all come from the {doc}`Kalman filter `. + +The state is the history of $Y_t$ that a finite-order VAR requires, and the +observation is the subvector $y_t$. + +Because $y_t$ is a subvector of the state and not a noisy signal about it, this +is a state space system with *no* measurement error. + +The main results are + +- $y_t$ obeys an **infinite-order** VAR whose coefficients we compute exactly, +- the innovation $a_t$ of the small system is a **one-sided distributed lag** of + the innovations $\varepsilon_t$ of the large system, +- that distributed lag has a wide coefficient matrix at every lag, so + $\varepsilon_t$ *cannot* be recovered from the history of $y_t$, +- the forecast error variance of the small system exceeds that of the large one + by a quantity we compute. + +Two special cases where the small VAR stays finite-order are identified and verified. + +This lecture generalizes the worked example that used to close +{doc}`kalman_filter_var`. + +Let's start with imports. + +```{code-cell} ipython3 +import matplotlib.pyplot as plt +import numpy as np +import quantecon as qe + +plt.rcParams['figure.figsize'] = (10, 5) +np.set_printoptions(precision=4, suppress=True) +``` + +## The large system + +Let $Y_t$ be $n \times 1$ and suppose it obeys the $m$th order vector autoregression + +$$ +Y_{t+1} = A_1 Y_t + A_2 Y_{t-1} + \cdots + A_m Y_{t-m+1} + \varepsilon_{t+1}, +\qquad +\mathbb{E}\, \varepsilon_t \varepsilon_t^\top = V , +$$ (eq:vs_bigvar) + +where $\{\varepsilon_t\}$ is a serially uncorrelated sequence with +$\mathbb{E}[\varepsilon_{t+1} \mid Y_t, Y_{t-1}, \ldots] = 0$. + +Stack $m$ lags into the $nm \times 1$ state vector + +$$ +X_t = \begin{pmatrix} Y_t \\ Y_{t-1} \\ \vdots \\ Y_{t-m+1}\end{pmatrix} . +$$ + +Then {eq}`eq:vs_bigvar` becomes the first-order **companion form** + +$$ +X_{t+1} = A X_t + C \varepsilon_{t+1}, +\qquad +A = \begin{pmatrix} +A_1 & A_2 & \cdots & A_{m-1} & A_m \\ +I & 0 & \cdots & 0 & 0 \\ +0 & I & \cdots & 0 & 0 \\ +\vdots & & \ddots & & \vdots \\ +0 & 0 & \cdots & I & 0 +\end{pmatrix}, +\qquad +C = \begin{pmatrix} I_n \\ 0 \\ \vdots \\ 0 \end{pmatrix} . +$$ (eq:vs_companion) + +We write $A$ without a subscript for the companion matrix and $A_1, \ldots, A_m$ +with subscripts for the VAR coefficient matrices. + +Let $J = \begin{pmatrix} I_n & 0 & \cdots & 0\end{pmatrix}$ be the $n \times nm$ +matrix that reads $Y_t$ off the state, so that $Y_t = J X_t$ and $C = J^\top$. + +We assume all eigenvalues of $A$ are strictly inside the unit circle, so +$\{Y_t\}$ is covariance stationary. + +## The small system + +The econometrician observes + +$$ +y_t = S_y Y_t = G X_t, +\qquad +G = S_y J , +$$ (eq:vs_obs) + +where $S_y$ is $n_y \times n$ with $n_y < n$. + +Usually $S_y$ picks out $n_y$ coordinates of $Y_t$, but nothing below requires +that, so linear combinations are allowed too. + +Equations {eq}`eq:vs_companion` and {eq}`eq:vs_obs` are a state space system of the +form studied in {doc}`kalman_filter_var`, with shock loading $C$, observation +matrix $G$, and measurement error covariance + +$$ +R = 0 . +$$ + +The observation is *exact*, but it is a strict subvector of the state, so the +econometrician still faces a filtering problem. + +Let $\Sigma$ solve the steady-state Riccati equation + +$$ +\Sigma = A \Sigma A^\top + C V C^\top + - A \Sigma G^\top \bigl(G \Sigma G^\top\bigr)^{-1} G \Sigma A^\top +$$ (eq:vs_riccati) + +with associated Kalman gain + +$$ +K = A \Sigma G^\top \bigl(G \Sigma G^\top\bigr)^{-1} . +$$ (eq:vs_gain) + +Here $\Sigma$ is the covariance matrix of $X_t - \mathbb{E}[X_t \mid y^{t-1}]$. + +```{note} +Because $X_t$ contains lags of $Y_t$ whose $S_y$ components the econometrician +has already seen exactly, $\Sigma$ is singular. + +That is harmless. + +What the Kalman gain {eq}`eq:vs_gain` requires is that the *innovation* +covariance $G \Sigma G^\top$ be nonsingular, which holds whenever no linear +combination of $y_t$ is perfectly predictable from $y^{t-1}$. +``` + +The innovation in the small system is + +$$ +a_t = y_t - \mathbb{E}[y_t \mid y^{t-1}] = G\bigl(X_t - \hat X_t\bigr), +\qquad +\Omega \equiv \mathbb{E}\, a_t a_t^\top = G \Sigma G^\top . +$$ (eq:vs_innov) + +## Four representations + +### The Wold representation + +The steady-state innovations representation derived in {doc}`kalman_filter_var` is + +$$ +\hat X_{t+1} = A \hat X_t + K a_t, +\qquad +y_t = G \hat X_t + a_t . +$$ (eq:vs_innovrep) + +Solving {eq}`eq:vs_innovrep` forward gives the moving average representation of +$y_t$ in terms of its own innovations, + +$$ +y_t = \sum_{h=0}^{\infty} \Psi_h\, a_{t-h}, +\qquad +\Psi_0 = I_{n_y}, +\qquad +\Psi_h = G A^{h-1} K \quad (h \geq 1) . +$$ (eq:vs_wold) + +### The vector autoregression + +Solving {eq}`eq:vs_innovrep` backward instead gives + +$$ +y_t = \sum_{j=1}^{\infty} B_j\, y_{t-j} + a_t, +\qquad +B_j = G (A - KG)^{j-1} K . +$$ (eq:vs_var) + +This is an infinite-order VAR, convergent because the eigenvalues of $A - KG$ +lie inside the unit circle. + +### Innovations of the small system in terms of the large one + +This is the question that motivates the lecture. + +Let $e_t = X_t - \hat X_t$ be the filtering error. + +Subtracting the Kalman recursion $\hat X_{t+1} = A \hat X_t + K a_t$ from the +state equation {eq}`eq:vs_companion` and using $a_t = G e_t$ gives + +$$ +e_{t+1} = (A - KG) e_t + C \varepsilon_{t+1} . +$$ (eq:vs_error) + +Solving {eq}`eq:vs_error` backward and premultiplying by $G$ yields the answer. + +```{prf:proposition} +:label: vs_prop_innov + +The innovations of the small system are the one-sided distributed lag + +$$ +a_t = \sum_{j=0}^{\infty} \Gamma_j\, \varepsilon_{t-j}, +\qquad +\Gamma_j = G (A - KG)^j C , +$$ (eq:vs_gamma) + +of the innovations of the large system, with leading coefficient + +$$ +\Gamma_0 = G C = S_y J J^\top = S_y . +$$ +``` + +Since each $\Gamma_j$ is $n_y \times n$ with $n_y < n$, the map from +$\{\varepsilon_t\}$ to $\{a_t\}$ has no inverse. + +Knowing the entire history of $y_t$ is not enough to recover $\varepsilon_t$. + +### The structural moving average + +For comparison, iterating {eq}`eq:vs_companion` gives $y_t$ directly in terms of +the large system's innovations, + +$$ +y_t = \sum_{j=0}^{\infty} \Phi_j\, \varepsilon_{t-j}, +\qquad +\Phi_j = G A^j C . +$$ (eq:vs_phi) + +### A forecast error that contains past shocks + +Separating the $j = 0$ term in {eq}`eq:vs_gamma` and using $\Gamma_0 = S_y$ gives + +$$ +a_t = \underbrace{S_y \varepsilon_t}_{\text{full-information forecast error}} + \; + \; + \underbrace{\sum_{j=1}^{\infty} \Gamma_j\, \varepsilon_{t-j}}_{\text{shocks realized before } t} . +$$ (eq:vs_split) + +The first term is what it appears to be. + +Because $\mathbb{E}[\varepsilon_t \mid Y^{t-1}] = 0$, we have +$y_t - \mathbb{E}[y_t \mid Y^{t-1}] = S_y \varepsilon_t$, so $S_y \varepsilon_t$ +is the error made in forecasting $y_t$ by someone who observes the *entire* +history of $Y$. + +The second term deserves a pause, because at first sight it looks impossible. + +By construction $a_t$ is a forecast error, orthogonal to everything known at +$t-1$. + +Yet {eq}`eq:vs_split` says that $a_t$ loads on $\varepsilon_{t-1}, +\varepsilon_{t-2}, \ldots$, shocks that had already been realized by then. + +The two facts are consistent, and both are true: + +$$ +\mathbb{E}\, a_t\, \varepsilon_{t-j}^\top = \Gamma_j V \neq 0 +\quad (j \geq 1), +\qquad \text{while} \qquad +\mathbb{E}\, a_t\, y_{t-k}^\top = 0 +\quad (k \geq 1) . +$$ (eq:vs_orth) + +The resolution is that past shocks are known to the *large* system's +econometrician, not to the small one. + +The small econometrician's information set is $H(y^{t-1})$, the closed linear +span of $y_{t-1}, y_{t-2}, \ldots$, and $\varepsilon_{t-j}$ does not lie in it. + +Let $P_{t-1}$ denote projection onto $H(y^{t-1})$. + +Applying $P_{t-1}$ to {eq}`eq:vs_split`, using $P_{t-1} a_t = 0$ and +$P_{t-1}\varepsilon_t = 0$, delivers the identity + +$$ +\sum_{j=1}^{\infty} \Gamma_j\, P_{t-1}\varepsilon_{t-j} = 0 . +$$ (eq:vs_pred) + +So the distributed lag in {eq}`eq:vs_split` loads only on the parts of past +shocks that the small econometrician has *not yet* learned, + +$$ +a_t = S_y \varepsilon_t + + \sum_{j=1}^{\infty} \Gamma_j + \bigl(\varepsilon_{t-j} - P_{t-1}\varepsilon_{t-j}\bigr) . +$$ (eq:vs_unlearned) + +A shock that happened three quarters ago can still be news today, if the only +series you watch has not finished revealing it. + +That is the whole content of the filtering problem, and it is why $\Omega$ +exceeds $S_y V S_y^\top$. + +### What the coefficients $\Gamma_j$ are + +Substituting the moving average {eq}`eq:vs_phi` into the autoregression +{eq}`eq:vs_var` and matching the coefficient on $\varepsilon_{t-k}$ gives a +second formula for the same objects, + +$$ +\Gamma_k = \Phi_k - \sum_{j=1}^{k} B_j\, \Phi_{k-j} . +$$ (eq:vs_gamma_alt) + +So $\Gamma_k$ measures the failure of the small system's own autoregression to +reproduce the large system's $k$-lag response. + +The case $k = 1$ is worth writing out. + +Since $\Phi_0 = S_y$ and $\Phi_1 = G A C = S_y A_1$, + +$$ +\Gamma_1 = S_y A_1 - B_1 S_y . +$$ (eq:vs_gamma1) + +Suppose $S_y$ selects coordinates, and partition +$Y_t = (y_t^\top, \tilde y_t^\top)^\top$ as before. + +Then $B_1 S_y$ has zeros in the columns belonging to the dropped variables, so +{eq}`eq:vs_gamma1` reads + +$$ +\Gamma_1 = \begin{pmatrix} A_1^{yy} - B_1 & A_1^{y \tilde y}\end{pmatrix} . +$$ (eq:vs_gamma1_block) + +The loading of $a_t$ on last period's *omitted* shocks is exactly +$A_1^{y\tilde y}$, the block through which the omitted variables enter the +retained equations. + +Equation {eq}`eq:vs_gamma1_block` also previews +{prf:ref}`vs_prop_blockexog`: block exogeneity sets $A_1^{y\tilde y} = 0$, which +forces $B_1 = A_1^{yy}$ and hence $\Gamma_1 = 0$. + +### A factorization identity + +Representations {eq}`eq:vs_wold`, {eq}`eq:vs_gamma`, and {eq}`eq:vs_phi` describe +the same process, so with $\Psi(z) = \sum_h \Psi_h z^h$ and similarly for +$\Gamma$ and $\Phi$, + +$$ +\Phi(z) = \Psi(z)\, \Gamma(z) . +$$ (eq:vs_factor) + +Equation {eq}`eq:vs_factor` says the structural moving average operator factors +into the Wold operator of the small system times the innovation map. + +It is a sharp numerical check on everything above, and we use it as one. + +Matching variances at each lag also gives + +$$ +\Omega = \sum_{j=0}^{\infty} \Gamma_j V \Gamma_j^\top + = S_y V S_y^\top + \sum_{j=1}^{\infty} \Gamma_j V \Gamma_j^\top . +$$ (eq:vs_varloss) + +Every term in the second sum is positive semidefinite, so + +$$ +\Omega \succeq S_y V S_y^\top . +$$ + +The small system's one-step forecast error variance is never smaller than the +corresponding block of the large system's, and {eq}`eq:vs_varloss` says exactly +how much is lost. + +## Two cases where nothing is lost + +```{prf:proposition} +:label: vs_prop_full + +If $S_y = I_n$, then $\Gamma_0 = I_n$ and $\Gamma_j = 0$ for $j \geq 1$, so +$a_t = \varepsilon_t$, and {eq}`eq:vs_var` collapses to the original $m$th order +VAR {eq}`eq:vs_bigvar`. +``` + +Nothing is hidden, so the Wold innovations *are* the structural innovations. + +The second case is more interesting. + +Partition $Y_t = (y_t^\top, \tilde y_t^\top)^\top$ and correspondingly + +$$ +A_k = \begin{pmatrix} A_k^{yy} & A_k^{y\tilde y} \\ + A_k^{\tilde y y} & A_k^{\tilde y \tilde y}\end{pmatrix}, +\qquad k = 1, \ldots, m . +$$ + +```{prf:proposition} +:label: vs_prop_blockexog + +Suppose $y_t$ is **block exogenous**, meaning $A_k^{y\tilde y} = 0$ for all $k$, +so that no lag of the omitted variables appears in the equations for $y$. + +Then $\Gamma_j = 0$ for $j \geq 1$, $a_t = S_y \varepsilon_t$, +$\Omega = S_y V S_y^\top$, and $y_t$ obeys the finite $m$th order VAR + +$$ +y_t = \sum_{k=1}^{m} A_k^{yy}\, y_{t-k} + a_t . +$$ +``` + +```{prf:proof} +Block exogeneity makes the $y$ rows of {eq}`eq:vs_bigvar` read +$y_{t+1} = \sum_k A_k^{yy} y_{t+1-k} + S_y \varepsilon_{t+1}$. + +The right side involves only lags of $y$, and $S_y \varepsilon_{t+1}$ is +orthogonal to the whole history $Y^t$ and hence to $y^t$. + +So this *is* the projection of $y_{t+1}$ on $y^t$, which identifies it as the +Wold representation. +``` + +Note what {prf:ref}`vs_prop_blockexog` does *not* require: $V$ need not be block +diagonal. + +Contemporaneous correlation between $\varepsilon^y$ and $\varepsilon^{\tilde y}$ +is fine. + +What matters for whether an omitted variable damages a VAR is Granger causality, +not contemporaneous correlation. + +## Code + +The class below packages everything. + +```{code-cell} ipython3 +class VARSubsystem: + """ + A VAR for Y and the implied representations for a subvector y = S_y Y. + + Y[t+1] = A_1 Y[t] + ... + A_m Y[t-m+1] + eps[t+1], E eps eps' = V + y[t] = S_y Y[t] + + Parameters + ---------- + A_list : list of (n, n) arrays, the VAR coefficient matrices A_1, ..., A_m + V : (n, n) array, covariance matrix of eps + S_y : (n_y, n) selector matrix + """ + + def __init__(self, A_list, V, S_y): + self.A_list = [np.atleast_2d(np.asarray(a, dtype=float)) for a in A_list] + self.V = np.atleast_2d(np.asarray(V, dtype=float)) + self.S_y = np.atleast_2d(np.asarray(S_y, dtype=float)) + n, m = self.A_list[0].shape[0], len(self.A_list) + self.n, self.m, self.n_y = n, m, self.S_y.shape[0] + + # companion form + self.A = np.zeros((n * m, n * m)) + self.A[:n] = np.hstack(self.A_list) + if m > 1: + self.A[n:, :n * (m - 1)] = np.eye(n * (m - 1)) + self.J = np.zeros((n, n * m)) + self.J[:, :n] = np.eye(n) + self.C = self.J.T + self.G = self.S_y @ self.J + self.Q = self.C @ self.V @ self.C.T + self._Sigma = self._K = None + + def companion_eigenvalues(self): + return np.linalg.eigvals(self.A) + + def stationary_filter(self): + """Steady-state (Sigma, K) from the Riccati equation with R = 0.""" + if self._Sigma is None: + A, G = self.A, self.G + R = np.zeros((self.n_y, self.n_y)) + Sigma = qe.solve_discrete_riccati(A.T, G.T, self.Q, R) + Omega = G @ Sigma @ G.T + self._Sigma = Sigma + self._K = A @ Sigma @ G.T @ np.linalg.inv(Omega) + return self._Sigma, self._K + + def innovation_cov(self): + """Omega = E a a', the one-step forecast error covariance of y.""" + Sigma, _ = self.stationary_filter() + return self.G @ Sigma @ self.G.T + + def wold(self, h_max=20): + """Psi[h] in y[t] = sum_h Psi[h] a[t-h]; Psi[0] = I.""" + _, K = self.stationary_filter() + Psi = np.empty((h_max + 1, self.n_y, self.n_y)) + Psi[0], P = np.eye(self.n_y), np.eye(self.A.shape[0]) + for h in range(1, h_max + 1): + Psi[h] = self.G @ P @ K + P = P @ self.A + return Psi + + def var_coefficients(self, h_max=20): + """B[j-1] in y[t] = sum_j B[j] y[t-j] + a[t], for j = 1, ..., h_max.""" + _, K = self.stationary_filter() + M = self.A - K @ self.G + B, P = np.empty((h_max, self.n_y, self.n_y)), np.eye(self.A.shape[0]) + for j in range(h_max): + B[j] = self.G @ P @ K + P = P @ M + return B + + def innovation_map(self, h_max=20): + """Gamma[j] in a[t] = sum_j Gamma[j] eps[t-j].""" + _, K = self.stationary_filter() + M = self.A - K @ self.G + Gamma, P = np.empty((h_max + 1, self.n_y, self.n)), np.eye(self.A.shape[0]) + for j in range(h_max + 1): + Gamma[j] = self.G @ P @ self.C + P = P @ M + return Gamma + + def structural_ma(self, h_max=20): + """Phi[j] in y[t] = sum_j Phi[j] eps[t-j].""" + Phi, P = np.empty((h_max + 1, self.n_y, self.n)), np.eye(self.A.shape[0]) + for j in range(h_max + 1): + Phi[j] = self.G @ P @ self.C + P = P @ self.A + return Phi + + def simulate(self, T, seed=0, burn=200): + """Simulate Y and the innovations eps that generated it.""" + rng = np.random.default_rng(seed) + L = np.linalg.cholesky(self.V) + eps = rng.standard_normal((T + burn, self.n)) @ L.T + X, Y = np.zeros(self.n * self.m), np.zeros((T + burn, self.n)) + for t in range(T + burn): + X = self.A @ X + self.C @ eps[t] + Y[t] = self.J @ X + return Y[burn:], eps[burn:] + + def filter_innovations(self, y_path): + """Recover a[t] from observed y by running the steady-state filter.""" + _, K = self.stationary_filter() + x_hat = np.zeros(self.A.shape[0]) + a = np.empty((len(y_path), self.n_y)) + for t in range(len(y_path)): + a[t] = y_path[t] - self.G @ x_hat + x_hat = self.A @ x_hat + K @ a[t] + return a + + +def convolve(Psi, Gamma, h_max): + """(Psi * Gamma)[h] = sum_{k=0}^{h} Psi[k] Gamma[h-k].""" + out = np.zeros((h_max + 1, Psi.shape[1], Gamma.shape[2])) + for h in range(h_max + 1): + for k in range(h + 1): + out[h] += Psi[k] @ Gamma[h - k] + return out +``` + +A single routine collects the diagnostics we want to see for every example. + +```{code-cell} ipython3 +def report(model, h_max=30, label=''): + """Print the identities that every subsystem must satisfy.""" + Sigma, K = model.stationary_filter() + A, G, V = model.A, model.G, model.V + resid = Sigma - (A @ Sigma @ A.T + model.Q + - A @ Sigma @ G.T @ np.linalg.inv(G @ Sigma @ G.T) + @ G @ Sigma @ A.T) + Psi = model.wold(h_max) + Gamma = model.innovation_map(h_max) + Phi = model.structural_ma(h_max) + Omega, Vy = model.innovation_cov(), model.S_y @ V @ model.S_y.T + + print(f'--- {label} (n = {model.n}, m = {model.m}, n_y = {model.n_y})') + print(f' max |eig| of companion A {np.max(abs(model.companion_eigenvalues())):.6f}') + print(f' max |eig| of A - KG {np.max(abs(np.linalg.eigvals(A - K @ G))):.6f}') + print(f' Riccati residual {np.abs(resid).max():.2e}') + print(f' |Gamma[0] - S_y| {np.abs(Gamma[0] - model.S_y).max():.2e}') + print(f' |Phi - Psi * Gamma| ' + f'{np.abs(convolve(Psi, Gamma, h_max) - Phi).max():.2e}') + Gamma_long = model.innovation_map(300) # the sum in (SS) is infinite + print(f' |Omega - sum Gamma V Gamma\'| ' + f'{np.abs(Omega - sum(Gamma_long[j] @ V @ Gamma_long[j].T for j in range(301))).max():.2e}') + print(f' max |Gamma[j]|, j >= 1 {np.abs(Gamma[1:]).max():.3e}') + print(f' det Omega / det S_y V S_y\' ' + f'{np.linalg.det(Omega) / np.linalg.det(Vy):.4f}') + return Psi, Gamma, Phi +``` + +## Example 1: a bivariate VAR(2) + +Two observable series $r_t$ and $z_t$ obey the VAR(2) + +$$ +\begin{pmatrix} r_{t+1} \\ z_{t+1}\end{pmatrix} += A_1 \begin{pmatrix} r_t \\ z_t \end{pmatrix} ++ A_2 \begin{pmatrix} r_{t-1} \\ z_{t-1}\end{pmatrix} ++ \varepsilon_{t+1}, +$$ + +with + +$$ +A_1 = \begin{pmatrix} 0.80 & 0.75 \\ 0 & 0.75 \end{pmatrix}, +\qquad +A_2 = \begin{pmatrix} 0.05 & -0.72 \\ 0 & 0.20 \end{pmatrix}, +\qquad +V = I_2 . +$$ + +Note that $z$ is block exogenous: its equation contains no lags of $r$. + +But $r$ is *not*, since $z$ enters the $r$ equation with both lags. + +So dropping $z$ should matter, while dropping $r$ should not. + +```{code-cell} ipython3 +A1 = np.array([[0.80, 0.75], + [0.00, 0.75]]) +A2 = np.array([[0.05, -0.72], + [0.00, 0.20]]) +V2 = np.eye(2) + +S_both = np.eye(2) # observe (r, z) +S_r = np.array([[1.0, 0.0]]) # observe r only +S_z = np.array([[0.0, 1.0]]) # observe z only + +mod_both = VARSubsystem([A1, A2], V2, S_both) +mod_r = VARSubsystem([A1, A2], V2, S_r) +mod_z = VARSubsystem([A1, A2], V2, S_z) + +Psi_both, Gam_both, _ = report(mod_both, label='observe r and z') +print() +Psi_r, Gam_r, _ = report(mod_r, label='observe r only') +print() +Psi_z, Gam_z, _ = report(mod_z, label='observe z only') +``` + +Every identity holds to machine precision. + +The three cases differ exactly as the propositions predict. + +Observing both variables gives $\Gamma_j = 0$ for $j \geq 1$, so +$a_t = \varepsilon_t$, as {prf:ref}`vs_prop_full` requires. + +Observing only $z$, which is block exogenous, also gives $\Gamma_j = 0$ for +$j \geq 1$, so $a_t = \varepsilon_{z,t}$, as {prf:ref}`vs_prop_blockexog` +requires. + +Observing only $r$ is different. + +Here $\Gamma_j \neq 0$ for $j \geq 1$, and the ratio of forecast error variances +reports how much the $r$-only econometrician loses. + +```{code-cell} ipython3 +print('observe r only:') +print(f' Omega = {mod_r.innovation_cov()[0, 0]:.4f}') +print(f' S_y V S_y\' = {(S_r @ V2 @ S_r.T)[0, 0]:.4f}') +print('\n Gamma[j] for j = 0, ..., 5 (rows: response of a to eps_r, eps_z)') +print(Gam_r[:6, 0, :]) +``` + +The forecast error variance is over 50 percent larger than $V_{11}$. + +The extra variance is entirely attributable to past $\varepsilon_z$ shocks that +the econometrician sees only through their effect on $r$. + +### The VAR for the subsystem + +```{code-cell} ipython3 +--- +mystnb: + figure: + caption: VAR coefficients for the subsystem + name: fig-vs-varcoef +--- +B_r = mod_r.var_coefficients(12) + +fig, ax = plt.subplots() +ax.stem(np.arange(1, 13), B_r[:, 0, 0], basefmt=' ') +ax.axhline(0, color='k', lw=0.6) +ax.set_xlabel('lag $j$') +ax.set_ylabel('$B_j$') +ax.set_title(r'Population VAR coefficients for $r_t$ when only $r$ is observed') +fig.tight_layout() +plt.show() + +print('B_1, ..., B_6:', np.round(B_r[:6, 0, 0], 5)) +print('\nfor comparison, A_1[0,0] and A_2[0,0]:', A1[0, 0], A2[0, 0]) +``` + +The infinite-order VAR for $r$ alone is dominated by two lags, but neither +coefficient equals the corresponding coefficient in the bivariate system. + +Dropping $z$ does not simply delete the $z$ columns of the VAR; it changes what +is left. + +### Wold impulse responses + +```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Wold responses when both variables are observed + name: fig-vs-wold-both +--- +H = 25 +h = np.arange(H + 1) +Psi_both = mod_both.wold(H) + +fig, axes = plt.subplots(2, 2, figsize=(10, 6), sharex=True) +names, shocks = [r'$r_t$', r'$z_t$'], [r'$a_{r,t}$', r'$a_{z,t}$'] +for i in range(2): + for j in range(2): + axes[i, j].plot(h, Psi_both[:, i, j], lw=2) + axes[i, j].axhline(0, color='k', lw=0.6, ls='--') + axes[i, j].set_title(f'{names[i]} to {shocks[j]}', fontsize=10) + if i == 1: + axes[i, j].set_xlabel('horizon $h$') +fig.suptitle('Wold responses, both variables observed') +fig.tight_layout() +plt.show() +``` + +Because $a_t = \varepsilon_t$ here, these Wold responses coincide with the +structural responses of the bivariate VAR. + +```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Wold response when only r is observed + name: fig-vs-wold-r +--- +Psi_r = mod_r.wold(H) +Phi_r = mod_r.structural_ma(H) + +fig, ax = plt.subplots() +ax.plot(h, Psi_r[:, 0, 0], lw=2, label=r'$\Psi_h$: $r_t$ to its own innovation $a_t$') +ax.plot(h, Phi_r[:, 0, 0], lw=2, ls='--', + label=r'$\Phi_h$: $r_t$ to $\varepsilon_{r,t}$') +ax.plot(h, Phi_r[:, 0, 1], lw=2, ls=':', + label=r'$\Phi_h$: $r_t$ to $\varepsilon_{z,t}$') +ax.axhline(0, color='k', lw=0.6) +ax.set_xlabel('horizon $h$') +ax.set_ylabel('response') +ax.set_title(r'Wold versus structural responses of $r_t$') +ax.legend() +fig.tight_layout() +plt.show() +``` + +The Wold response to $a_t$ is not the response to either structural shock. + +It is a blend, and {eq}`eq:vs_factor` says exactly how the blending works. + +### Checking the innovation map by simulation + +Proposition {prf:ref}`vs_prop_innov` is a statement about population objects. + +We check it on a long simulated sample by comparing the innovations that the +Kalman filter extracts from the observed $r$ series with the distributed lag +$\sum_j \Gamma_j \varepsilon_{t-j}$ built from the shocks that generated it. + +```{code-cell} ipython3 +T = 100_000 +Y_sim, eps_sim = mod_r.simulate(T, seed=1) +y_sim = Y_sim @ S_r.T + +a_filtered = mod_r.filter_innovations(y_sim) + +J_lag = 400 # the distributed lag is infinite, so truncate generously +Gam_long = mod_r.innovation_map(J_lag) +a_theory = np.zeros_like(a_filtered) +for j in range(J_lag + 1): + a_theory[j:] += eps_sim[:T - j] @ Gam_long[j].T + +burn = J_lag + 1 +print(f'correlation of the two series ' + f'{np.corrcoef(a_filtered[burn:, 0], a_theory[burn:, 0])[0, 1]:.8f}') +print(f'max absolute difference ' + f'{np.abs(a_filtered[burn:] - a_theory[burn:]).max():.2e}') +print(f'sample variance of a {a_filtered[burn:, 0].var():.4f}') +print(f'population Omega {mod_r.innovation_cov()[0, 0]:.4f}') +print(f'sample corr(a_t, a_(t-1)) ' + f'{np.corrcoef(a_filtered[burn + 1:, 0], a_filtered[burn:-1, 0])[0, 1]:.4f}') +``` + +The two constructions of $a_t$ agree to the accuracy of the truncated +distributed lag, the sample variance of $a_t$ matches $\Omega$, and $a_t$ is +serially uncorrelated. + +### The two orthogonality facts + +Now we check {eq}`eq:vs_orth` directly, by running two regressions on the +simulated data. + +The first regresses $a_t$ on current and lagged *structural* shocks, which the +small econometrician cannot see. + +The second regresses $a_t$ on lagged *observations*, which are all that the small +econometrician can see. + +```{code-cell} ipython3 +P_lags = 4 +Z_eps = np.column_stack([eps_sim[P_lags - l:T - l] for l in range(P_lags + 1)]) +target = a_filtered[P_lags:, 0] +b_eps = np.linalg.lstsq(Z_eps, target, rcond=None)[0].reshape(P_lags + 1, 2) +fit = Z_eps @ b_eps.ravel() +r2_eps = 1 - ((target - fit) ** 2).sum() / target.var() / len(target) + +print('OLS of a_t on eps_t, ..., eps_{t-4}:') +print(b_eps) +print('population Gamma_0, ..., Gamma_4:') +print(Gam_long[:P_lags + 1, 0, :]) +print(f'max discrepancy {np.abs(b_eps - Gam_long[:P_lags + 1, 0, :]).max():.2e}') +print(f'R^2 = {r2_eps:.6f}') +``` + +```{code-cell} ipython3 +Q_lags = 8 +Z_y = np.column_stack([y_sim[Q_lags - l - 1:T - l - 1, 0] for l in range(Q_lags)]) +tgt = a_filtered[Q_lags:, 0] +b_y = np.linalg.lstsq(Z_y, tgt, rcond=None)[0] +resid = tgt - Z_y @ b_y +r2_y = 1 - (resid ** 2).sum() / tgt.var() / len(tgt) + +print('OLS of a_t on y_{t-1}, ..., y_{t-8}:') +print(np.round(b_y, 5)) +print(f'R^2 = {r2_y:.6f}') +``` + +The first regression recovers the $\Gamma_j$ and fits almost perfectly, the +small shortfall coming only from truncating the distributed lag at four lags. + +So $a_t$ really is built out of shocks stretching back before $t$. + +The second explains essentially nothing, confirming that $a_t$ is nonetheless +orthogonal to the small econometrician's information set. + +We can see why by asking how much of each past shock the small econometrician has +managed to learn. + +```{code-cell} ipython3 +print('R^2 from projecting a structural shock on y_{t-1}, ..., y_{t-8}') +for lag in [0, 1, 2, 3]: + r2s = [] + for k in range(2): + shock = eps_sim[Q_lags - lag:T - lag, k] + c = np.linalg.lstsq(Z_y, shock, rcond=None)[0] + e = shock - Z_y @ c + r2s.append(1 - (e ** 2).sum() / shock.var() / len(shock)) + print(f' eps_(t-{lag}): eps_r {r2s[0]:6.4f} eps_z {r2s[1]:6.4f}') +``` + +The current shock $\varepsilon_t$ is entirely unpredictable from $y^{t-1}$, as it +must be. + +Shocks from two or more periods back are substantially learned. + +The interesting row is $\varepsilon_{t-1}$: its $r$ component is largely known, +while its $z$ component is *completely* unknown, because $z_{t-1}$ reaches $r$ +only with a one-period lag and so has not yet shown up anywhere in $y^{t-1}$. + +That is why {eq}`eq:vs_gamma1_block` gives $a_t$ a loading on +$\varepsilon_{z,t-1}$ equal to the full structural coefficient +$A_1^{y\tilde y} = 0.75$, while its loading on $\varepsilon_{r,t-1}$ is only the +much smaller residual $A_1^{yy} - B_1$. + +```{code-cell} ipython3 +B1 = mod_r.var_coefficients(1)[0] +print(f'Gamma_1 = {Gam_long[1, 0, :]}') +print(f'[A1[0,0] - B_1, A1[0,1]] = ' + f'[{A1[0, 0] - B1[0, 0]:.4f}, {A1[0, 1]:.4f}]') + +Phi_r = mod_r.structural_ma(6) +B_r6 = mod_r.var_coefficients(6) +recursion = np.array([Phi_r[k] - sum(B_r6[j - 1] @ Phi_r[k - j] + for j in range(1, k + 1)) + for k in range(7)]) +print(f'\nmax |Gamma_k - (Phi_k - sum_j B_j Phi_(k-j))| = ' + f'{np.abs(recursion - Gam_long[:7]).max():.2e}') +``` + +## Example 2: an omitted interest rate + +Now a trivariate VAR(1) in output growth $g_t$, inflation $\pi_t$, and an +interest rate $i_t$, from which the econometrician drops $i_t$. + +We contrast two coefficient matrices that differ *only* in whether the interest +rate feeds back onto $g$ and $\pi$. + +$$ +A_1^{\text{exog}} = +\begin{pmatrix} +0.60 & 0.10 & 0.00 \\ +0.15 & 0.55 & 0.00 \\ +0.30 & 0.40 & 0.70 +\end{pmatrix}, +\qquad +A_1^{\text{fb}} = +\begin{pmatrix} +0.60 & 0.10 & -0.35 \\ +0.15 & 0.55 & 0.25 \\ +0.30 & 0.40 & 0.70 +\end{pmatrix} . +$$ + +The shock covariance matrix $V$ is the same in both, and it is *not* diagonal, so +the interest rate innovation is contemporaneously correlated with the other two. + +```{code-cell} ipython3 +V3 = np.array([[0.36, 0.05, 0.02], + [0.05, 0.25, 0.06], + [0.02, 0.06, 0.16]]) +S_gpi = np.array([[1.0, 0.0, 0.0], + [0.0, 1.0, 0.0]]) + +A_exog = np.array([[0.60, 0.10, 0.00], + [0.15, 0.55, 0.00], + [0.30, 0.40, 0.70]]) +A_fb = np.array([[0.60, 0.10, -0.35], + [0.15, 0.55, 0.25], + [0.30, 0.40, 0.70]]) + +mod_exog = VARSubsystem([A_exog], V3, S_gpi) +mod_fb = VARSubsystem([A_fb], V3, S_gpi) + +_, Gam_exog, _ = report(mod_exog, label='i is block exogenous') +print() +_, Gam_fb, _ = report(mod_fb, label='i feeds back') +``` + +The block exogenous case behaves exactly as {prf:ref}`vs_prop_blockexog` says it +must, despite the correlated shocks. + +The feedback case does not. + +```{code-cell} ipython3 +print('block exogenous: B_1 versus the (g, pi) block of A_1') +print(mod_exog.var_coefficients(3)[0]) +print(A_exog[:2, :2]) +print(f' max |B_j| for j >= 2: {np.abs(mod_exog.var_coefficients(12)[1:]).max():.2e}') + +print('\nfeedback: B_1 versus the (g, pi) block of A_1') +print(mod_fb.var_coefficients(3)[0]) +print(A_fb[:2, :2]) +print(f' max |B_j| for j >= 2: {np.abs(mod_fb.var_coefficients(12)[1:]).max():.4f}') +``` + +With block exogeneity the subsystem VAR is *exactly* the corresponding block of +the large VAR, and it stops at one lag. + +With feedback the one-lag coefficients are distorted and higher-order terms +appear. + +The effect on the coefficient of lagged inflation in the output growth equation +is worth noticing: a positive number in the large system becomes a negative one +in the subsystem. + +```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Innovation map coefficients with and without feedback + name: fig-vs-gamma +--- +J_max = 10 +labels = [r'$\varepsilon_g$', r'$\varepsilon_\pi$', r'$\varepsilon_i$'] +fig, axes = plt.subplots(2, 2, figsize=(11, 6), sharex=True) +for col, (Gam, ttl) in enumerate([(Gam_exog, 'block exogenous'), + (Gam_fb, 'feedback')]): + for row, obs in enumerate([r'$a_g$', r'$a_\pi$']): + ax = axes[row, col] + for k in range(3): + ax.plot(np.arange(J_max + 1), Gam[:J_max + 1, row, k], + marker='o', ms=3, lw=1.5, label=labels[k]) + ax.axhline(0, color='k', lw=0.6) + ax.set_title(f'{obs}, {ttl}', fontsize=10) + if row == 1: + ax.set_xlabel('lag $j$') + if row == 0 and col == 0: + ax.legend(fontsize=8) +axes[0, 0].set_ylabel(r'$\Gamma_j$') +axes[1, 0].set_ylabel(r'$\Gamma_j$') +fig.suptitle(r'Coefficients $\Gamma_j$ in $a_t = \sum_j \Gamma_j \varepsilon_{t-j}$') +fig.tight_layout() +plt.show() +``` + +In the left column only the $j = 0$ coefficients are nonzero, and they equal the +rows of $S_y$. + +In the right column the omitted interest rate shock $\varepsilon_i$ leaks into +the observed innovations at every lag. + +## Summary + +A finite-order VAR for $Y_t$ implies, for any subvector $y_t = S_y Y_t$, an +infinite-order VAR whose coefficients $B_j = G(A - KG)^{j-1}K$ come from the +steady-state Kalman filter for the companion system. + +The innovations of the small system are a one-sided distributed lag +$a_t = \sum_j \Gamma_j \varepsilon_{t-j}$ of the innovations of the large system, +with $\Gamma_j = G(A - KG)^j C$ and $\Gamma_0 = S_y$. + +Because $\Gamma_j$ is wide, that map cannot be inverted, which is a precise +statement of the informational deficiency of a subsystem VAR. + +The price is measured by $\Omega - S_y V S_y^\top = \sum_{j \geq 1} \Gamma_j V \Gamma_j^\top$. + +The price is zero when everything is observed, and also when the retained block +is block exogenous, in which case the subsystem VAR is exactly the corresponding +block of the original one. + +## Exercises + +```{exercise-start} +:label: vs_ex1 +``` + +Take the trivariate system of Example 2 and put the feedback of the interest rate +onto $(g, \pi)$ under your control by writing + +$$ +A_1(\theta) = A_1^{\text{exog}} + \theta \begin{pmatrix} 0 & 0 & -0.35 \\ +0 & 0 & 0.25 \\ 0 & 0 & 0 \end{pmatrix} . +$$ + +For $\theta$ on a grid from $0$ to $1.5$, plot + +1. $\det \Omega / \det(S_y V S_y^\top)$, the information lost by dropping $i_t$, +2. $\max_{j \geq 1} |\Gamma_j|$, the size of the leakage of past shocks into the + observed innovations. + +Explain the shape you find at $\theta = 0$. + +```{exercise-end} +``` + +```{solution-start} vs_ex1 +:class: dropdown +``` + +Here is one solution: + +```{code-cell} ipython3 +E = np.zeros((3, 3)) +E[0, 2], E[1, 2] = -0.35, 0.25 + +thetas = np.linspace(0, 1.5, 31) +det_ratio, leak = [], [] +for th in thetas: + mod = VARSubsystem([A_exog + th * E], V3, S_gpi) + Om = mod.innovation_cov() + det_ratio.append(np.linalg.det(Om) + / np.linalg.det(S_gpi @ V3 @ S_gpi.T)) + leak.append(np.abs(mod.innovation_map(40)[1:]).max()) + +fig, axes = plt.subplots(1, 2, figsize=(11, 4)) +axes[0].plot(thetas, det_ratio, lw=2) +axes[0].set(xlabel=r'$\theta$', ylabel='determinant ratio', + title='information lost by dropping $i_t$') +axes[1].plot(thetas, leak, lw=2, color='C1') +axes[1].set(xlabel=r'$\theta$', ylabel=r'$\max_{j \geq 1} |\Gamma_j|$', + title='leakage of past shocks into $a_t$') +for ax in axes: + ax.axhline(ax.get_ylim()[0], color='k', lw=0.6) +fig.tight_layout() +plt.show() + +print(f'at theta = 0: ratio = {det_ratio[0]:.6f}, leakage = {leak[0]:.2e}') +print(f'at theta = 1: ratio = {det_ratio[20]:.4f}, leakage = {leak[20]:.4f}') +``` + +Both curves start at their minima, exactly zero leakage and a determinant ratio +of exactly one. + +That is {prf:ref}`vs_prop_blockexog`: at $\theta = 0$ the interest rate does not +Granger cause $(g, \pi)$, so dropping it costs nothing even though its +innovation is contemporaneously correlated with the others. + +Both measures rise as the feedback strengthens. + +```{solution-end} +``` + +```{exercise-start} +:label: vs_ex2 +``` + +Return to the bivariate system, but replace $A_1$ and $A_2$ by the single matrix + +$$ +A_1 = \begin{pmatrix} 0.5 & 0.6 \\ 0 & \rho_z \end{pmatrix}, +\qquad V = I_2 , +$$ + +so that $\rho_z$ controls the persistence of the omitted variable $z$. + +Observing $r$ only, report for $\rho_z \in \{0.2, 0.5, 0.75, 0.9, 0.95, 0.99\}$ + +1. the forecast error variance $\Omega$, +2. the number of lags needed before $\sum_{j > p} |B_j| < 10^{-3}$. + +What does a persistent omitted variable do to the VAR that an econometrician +should fit? + +```{exercise-end} +``` + +```{solution-start} vs_ex2 +:class: dropdown +``` + +Here is one solution: + +```{code-cell} ipython3 +print(' rho_z Omega lags needed') +for rho_z in [0.2, 0.5, 0.75, 0.9, 0.95, 0.99]: + mod = VARSubsystem([np.array([[0.5, 0.6], [0.0, rho_z]])], + np.eye(2), np.array([[1.0, 0.0]])) + B = mod.var_coefficients(400) + tails = [np.abs(B[p:]).sum() for p in range(400)] + p_need = next(p for p in range(400) if tails[p] < 1e-3) + print(f' {rho_z:5.2f} {mod.innovation_cov()[0, 0]:6.4f} {p_need:3d}') +``` + +Both columns rise with $\rho_z$. + +A persistent omitted variable both inflates the forecast error variance and +lengthens the autoregression, because the econometrician must reach further back +to extract the same information about $z$ from the history of $r$. + +An empirical VAR with too few lags will therefore be worst exactly where the +missing variable is most persistent. + +```{solution-end} +``` + +```{exercise-start} +:label: vs_ex3 +``` + +The coefficients $B_j$ are population objects. + +Simulate $T = 200{,}000$ observations of the bivariate system of Example 1, +retain only $r_t$, and fit finite-order autoregressions of orders +$p = 1, 2, 4, 8$ by ordinary least squares. + +Compare the estimates with $B_1, \ldots, B_p$ and the residual variance with +$\Omega$. + +Which order is too short, and what does fitting too short a lag length do to the +first coefficient? + +```{exercise-end} +``` + +```{solution-start} vs_ex3 +:class: dropdown +``` + +Here is one solution: + +```{code-cell} ipython3 +Y_big, _ = mod_r.simulate(200_000, seed=3) +r_series = Y_big[:, 0] +B_pop = mod_r.var_coefficients(8)[:, 0, 0] + +for p in [1, 2, 4, 8]: + X = np.column_stack([r_series[p - 1 - l:len(r_series) - 1 - l] + for l in range(p)]) + zz = r_series[p:] + b_hat = np.linalg.lstsq(X, zz, rcond=None)[0] + resid = zz - X @ b_hat + print(f'p = {p}') + print(f' OLS {np.round(b_hat, 4)}') + print(f' population {np.round(B_pop[:p], 4)}') + print(f' residual variance {resid.var():.4f} Omega {mod_r.innovation_cov()[0, 0]:.4f}') +``` + +An AR(1) is too short. + +Its single coefficient is pulled well above $B_1$, because it has to stand in for +the omitted second lag, and its residual variance exceeds $\Omega$. + +From $p = 2$ onward the estimates track the population coefficients and the +residual variance settles on $\Omega$, which matches the finding above that +$B_j$ is negligible beyond the second lag in this example. + +```{solution-end} +``` From d1691647f4814f1b0092a64d379e2312dce14085 Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Sat, 1 Aug 2026 11:11:00 +1000 Subject: [PATCH 2/3] [sargent_surico] Update arviz calls for the 1.x API The lecture installs arviz unpinned, and `pip install arviz` now resolves 1.2.0 -- the rewrite that dispatches to arviz-base/arviz-stats/arviz-plots. Two calls still used the 0.x API. `az.from_dict` took the posterior group as its first positional argument in 0.x, so a flat {var_name: array} mapping worked. In 1.x the first argument is {group_name: {var_name: array}}, so each entry of FREE was read as a group and dict_to_dataset() was handed a raw ndarray, failing the HTML build with `AttributeError: 'numpy.ndarray' object has no attribute 'items'`. `az.summary` no longer emits hdi_3%/hdi_97%; it defaults to an 89% ETI. Ask for ci_kind='hdi', ci_prob=0.94 to keep the interval the 0.x default gave, and select the new hdi94_lb/hdi94_ub columns. Execution stopped at the from_dict cell above, so this second break never reached the CI log. Co-Authored-By: Claude Opus 5 (1M context) --- lectures/sargent_surico.md | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/lectures/sargent_surico.md b/lectures/sargent_surico.md index 205d84bd1..b9d706818 100644 --- a/lectures/sargent_surico.md +++ b/lectures/sargent_surico.md @@ -1284,7 +1284,8 @@ A correlated chain of length $N$ is worth fewer than $N$ independent draws, and ```{code-cell} ipython3 import arviz as az -rwmh_idata = az.from_dict({n: kept[:, j][None, :] for j, n in enumerate(FREE)}) +rwmh_idata = az.from_dict( + {'posterior': {n: kept[:, j][None, :] for j, n in enumerate(FREE)}}) rwmh_summary = az.summary(rwmh_idata, var_names=FREE) print(rwmh_summary[['mean', 'sd', 'ess_bulk']].to_string()) print(f'\n{len(kept)} retained draws out of {N_DRAWS}, in {rwmh_seconds:.0f} seconds') @@ -1637,9 +1638,10 @@ print(f'divergences ' ```{code-cell} ipython3 nuts_idata = az.from_numpyro(mcmc) -nuts_summary = az.summary(nuts_idata, var_names=FREE) +nuts_summary = az.summary(nuts_idata, var_names=FREE, + ci_kind='hdi', ci_prob=0.94) nuts_kept = np.column_stack([np.asarray(mcmc.get_samples()[n]) for n in FREE]) -nuts_summary[['mean', 'sd', 'hdi_3%', 'hdi_97%', 'ess_bulk', 'r_hat']] +nuts_summary[['mean', 'sd', 'hdi94_lb', 'hdi94_ub', 'ess_bulk', 'r_hat']] ``` ### A warning sign in the diagnostics From d812aa6137c0d6ed9cfebee61f63c585b3f3ba3a Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Sat, 1 Aug 2026 12:54:52 +1000 Subject: [PATCH 3/3] [sargent_surico] Pin jax to the CPU The NUTS cell exceeds the 2400s per-cell execution timeout set in lectures/_config.yml when it runs on the GPU CI runner, so the notebook dies with a CellTimeoutError and the HTML build fails. Executed locally on CPU, the same cell takes ~160s and the whole lecture takes ~5 minutes across all 52 code cells. Scaling by the ratio between local and CI timings for the non-jax part of the lecture puts the CPU cost of that cell at roughly 10 minutes on the runner, well inside the timeout. Why the GPU is so much slower here has not been confirmed. Pinning the platform is the change that makes the lecture build; the underlying cause is worth a separate look. Co-Authored-By: Claude Opus 5 (1M context) --- lectures/sargent_surico.md | 1 + 1 file changed, 1 insertion(+) diff --git a/lectures/sargent_surico.md b/lectures/sargent_surico.md index b9d706818..7b0534c13 100644 --- a/lectures/sargent_surico.md +++ b/lectures/sargent_surico.md @@ -1393,6 +1393,7 @@ from jax import lax from numpyro.infer import MCMC, NUTS jax.config.update('jax_enable_x64', True) +jax.config.update('jax_platform_name', 'cpu') U, NJ = 8, 9 # y = [pi, x, dm, R, e, a, chi, z, u]