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..7b0534c13 --- /dev/null +++ b/lectures/sargent_surico.md @@ -0,0 +1,2099 @@ +--- +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( + {'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') +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) +jax.config.update('jax_platform_name', 'cpu') + +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, + 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', 'hdi94_lb', 'hdi94_ub', '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} +```