From efc332431cfc2c6da560e638e8267ffda3fea9c7 Mon Sep 17 00:00:00 2001 From: thomassargent30 Date: Fri, 11 Sep 2026 15:56:47 -0400 Subject: [PATCH 1/2] [var_dmd, svd_intro, measurement_models] Correct math, add exercises var_dmd: reframe DMD as a rank-r VAR (principal components regression); replace the incorrect identity A = Phi Lambda Phi^+ with the exact modal representation A = Phi Lambda Psi, Psi = (W Lambda)^{-1} U_r^T; add a state-space form, long-run responses, a simulated example, and three exercises. Cite the published Sargent-Selvakumar-Yang (IER 2026) paper. svd_intro: fix the "not not" typo and other errors (four-subspaces code mislabeled V^T, Eckart-Young statement, reduced-SVD identities, PCA score dimensions and indices, polar decomposition); fix DecomAnalysis for tall matrices and use it; add two exercises bridging to var_dmd. measurement_models: use the innovation covariance Sigma_eta (not the unconditional R) in the quasi-differenced Kalman filter; fix the bracket in the consumption rule and shock timing; add an ordering-free Granger causality check, note the exact identity in filtered data, and add three exercises. Co-Authored-By: Claude Opus 5 (1M context) --- lectures/_static/quant-econ.bib | 49 ++ lectures/measurement_models.md | 409 +++++++++--- lectures/svd_intro.md | 649 +++++++++++++------ lectures/var_dmd.md | 1046 +++++++++++++++++-------------- 4 files changed, 1436 insertions(+), 717 deletions(-) diff --git a/lectures/_static/quant-econ.bib b/lectures/_static/quant-econ.bib index 1c082c164..86bf3c246 100644 --- a/lectures/_static/quant-econ.bib +++ b/lectures/_static/quant-econ.bib @@ -5356,3 +5356,52 @@ @article{ShakedSutton1982 pages = {3--13}, year = {1982} } + +@article{sargent2026dynamic, + title = {Dynamic Mode Decompositions and Vector Autoregressions}, + author = {Sargent, Thomas J. and Selvakumar, Yatheesan J. and Yang, Ziyue}, + journal = {International Economic Review}, + year = {2026}, + publisher = {Wiley Online Library} +} + +@techreport{SSY_CEX_2026, + author = {Sargent, Thomas J. and Selvakumar, Yatheesan J. and Yang, Ziyue}, + title = {Aggregate Shocks and Cross-Section Dynamics: Quantifying + Redistribution and Insurance in {US} Household Data}, + institution = {New York University}, + year = {2026} +} + +@article{anderson1951, + author = {Anderson, T. W.}, + title = {Estimating Linear Restrictions on Regression Coefficients for + Multivariate Normal Distributions}, + journal = {Annals of Mathematical Statistics}, + volume = {22}, + number = {3}, + pages = {327--351}, + year = {1951} +} + +@article{stock_watson2002, + author = {Stock, James H. and Watson, Mark W.}, + title = {Forecasting Using Principal Components from a Large Number of + Predictors}, + journal = {Journal of the American Statistical Association}, + volume = {97}, + number = {460}, + pages = {1167--1179}, + year = {2002} +} + +@article{Geweke1982, + author = {Geweke, John}, + title = {Measurement of Linear Dependence and Feedback Between Multiple + Time Series}, + journal = {Journal of the American Statistical Association}, + volume = {77}, + number = {378}, + pages = {304--313}, + year = {1982} +} \ No newline at end of file diff --git a/lectures/measurement_models.md b/lectures/measurement_models.md index ea4133a81..1746a53b5 100644 --- a/lectures/measurement_models.md +++ b/lectures/measurement_models.md @@ -41,9 +41,9 @@ If accurate observations on these time series are available, one can use that mapping to implement parameter estimation methods based either on the likelihood function or on the method of moments. -```{note} This is why econometrics estimation is often called an ''inverse'' problem, while -simulating a model for given parameter values is called a ''direct problem''. The direct problem -refers to the mapping we have just described, while the inverse problem involves somehow applying an ''inverse'' of that mapping to a data set that is treated as if it were one draw from the joint probability distribution described by the mapping. +```{note} This is why econometric estimation is often called an "inverse" problem, while +simulating a model for given parameter values is called a "direct" problem. The direct problem +refers to the mapping we have just described, while the inverse problem involves somehow applying an "inverse" of that mapping to a data set that is treated as if it were one draw from the joint probability distribution described by the mapping. ``` However, if only error-ridden data exist for the variables of interest, @@ -65,7 +65,7 @@ investment accelerator, the mechanism studied in these two quantecon lectures -- - In Model 1, the data collecting agency simply reports the error-ridden data that it collects. -- In Model 2, the data collection agents first collects error-ridden data that satisfy +- In Model 2, the data collection agency first collects error-ridden data that satisfy a classical errors-in-variables model, then filters the data, and reports the filtered objects. Although the two models have the same "deep parameters," they produce @@ -156,7 +156,7 @@ and $\theta_t$ is an endowment or technology shock following :label: shock_process a(L)\,\theta_t = \varepsilon_t, ``` -where $L$ is the backward shift (or 'lag') operator and $a(z) = 1 - a_1 z - a_2 z^2 - \cdots - a_r z^r$ having all its zeroes +where $L$ is the backward shift (or 'lag') operator, $\varepsilon_t$ is a scalar white noise, and $a(z) = 1 - a_1 z - a_2 z^2 - \cdots - a_r z^r$ has all its zeroes outside the unit circle. ### Optimal decision rule @@ -167,8 +167,10 @@ The optimal decision rule for $c_t$ is :label: opt_decision c_t = \frac{-\alpha}{f-1} + \left(1 - \frac{1}{\beta f^2}\right) + \left[ \frac{L - f^{-1} a(f^{-1})^{-1} a(L)}{L - f^{-1}}\,\theta_t - + f k_t, + + f k_t + \right], \qquad k_{t+1} = f k_t + \theta_t - c_t, ``` @@ -177,7 +179,7 @@ where $\alpha = u_1[1-(\beta f)^{-1}]/u_2$. Equations {eq}`shock_process` and {eq}`opt_decision` exhibit the cross-equation restrictions characteristic of rational expectations -models. +models: the parameters $a_1, \ldots, a_r$ of the shock process also appear in the decision rule for consumption. ### Net income and the accelerator @@ -211,13 +213,15 @@ measurement error. Assumption 2 is less important, affecting only various constants. -Under both assumptions, {eq}`opt_decision` simplifies to +Under both assumptions, $\alpha = 0$ and {eq}`opt_decision` simplifies to ```{math} :label: simple_crule -c_t = (1-f^{-1})\,\theta_t + (f-1)\,k_t. +c_t = (1-f^{-1})\,\theta_t + (f-1)\,k_t, ``` +so that the capital stock obeys $k_{t+1} = k_t + f^{-1} \theta_t$. + When {eq}`simple_crule`, {eq}`net_income`, and {eq}`tech_constraint` are combined, the optimal plan satisfies @@ -247,17 +251,21 @@ of income. This is the same mechanism that {cite:t}`Chow1968` documented empirically (see {doc}`chow_business_cycles`). -Equation {eq}`income_process` states that the first difference of disposable income is a +Equation {eq}`income_process` states that the first difference of net income, $y_{nt} - y_{n,t-1} = \theta_t - \beta \theta_{t-1}$, is a first-order moving average process with innovation equal to the innovation of the endowment shock $\theta_t$. As {cite:t}`Muth1960` showed, such a process is optimally forecast via a geometric distributed lag or "adaptive expectations" scheme. +Consumption itself is a martingale: {eq}`simple_crule` and $k_t - k_{t-1} = f^{-1}\theta_{t-1}$ imply $c_t - c_{t-1} = (1 - \beta)\theta_t$, the random walk property of consumption emphasized by {cite:t}`Hall1978`. + +So consumption, the capital stock, and net income all have unit roots, while investment $k_{t+1} - k_t = f^{-1}\theta_t$ is serially uncorrelated. + ### The accelerator puzzle When all variables are measured accurately and are driven by -the single shock $\theta_t$, the spectral density matrix of -$(c_t,\, k_{t+1}-k_t,\, y_{nt})$ has rank one at all frequencies. +the single shock $\theta_t$, the spectral density matrix of the covariance-stationary vector +$(c_t - c_{t-1},\, k_{t+1}-k_t,\, y_{nt} - y_{n,t-1})$ has rank one at all frequencies. Each variable is an invertible one-sided distributed lag of the same white noise, so no variable Granger-causes any other. @@ -288,10 +296,7 @@ series. Then + (v_{k,t+1} - v_{kt}), ``` -```{math} -:label: income_process_ma -y_{nt} = \theta_t + (1-\beta)(\theta_{t-1} + \theta_{t-2} + \cdots). -``` +where $y_{nt}$ is given by {eq}`income_process`. In this case income Granger-causes consumption and investment but is not Granger-caused by them. @@ -306,8 +311,7 @@ variances of the measurement errors. In this case, each observed series mixes the common signal $\theta_t$ with idiosyncratic measurement noise. -A series with lower measurement -error variance tracks $\theta_t$ more closely, so its innovations +A series whose measurement error is small **relative to the variation in the true series** tracks $\theta_t$ more closely, so its innovations contain more information about future values of the other series. Accordingly, in a forecast-error-variance decomposition, shocks to @@ -335,18 +339,20 @@ x_t = \begin{bmatrix} k_t \\ \theta_t \end{bmatrix}, z_t = \begin{bmatrix} y_{nt} \\ c_t \\ \Delta k_t \end{bmatrix}, ``` -so that the error-free data are described by the state-space system +where $\Delta k_t \equiv k_{t+1} - k_t$ denotes investment at time $t$. + +The error-free data are described by the state-space system ```{math} :label: true_ss \begin{aligned} -x_{t+1} &= A x_t + \varepsilon_t, \\ -z_t &= C x_t. +x_{t+1} &= A x_t + w_t, \\ +z_t &= C x_t, \end{aligned} ``` -where $\varepsilon_t = \begin{bmatrix} 0 \\ \theta_t \end{bmatrix}$ has -covariance $E \varepsilon_t \varepsilon_t^\top = Q$ and the matrices are +where $w_t = \begin{bmatrix} 0 \\ \theta_{t+1} \end{bmatrix}$ has +covariance $E w_t w_t^\top = Q$ and the matrices are ```{math} A = \begin{bmatrix} @@ -366,10 +372,14 @@ Q = \begin{bmatrix} \end{bmatrix}. ``` +Here we date the shock $w_t$ that arrives between $t$ and $t+1$ by $t$; this matches the timing convention for the measurement errors below. + $Q$ is singular because there is only one source of randomness $\theta_t$; the capital stock $k_t$ evolves deterministically given $\theta_t$. +Notice that the rows of $C$ satisfy $C_c + C_{\Delta k} = C_y$, which is the national income identity $c_t + \Delta k_t = y_{nt}$. + ```{code-cell} ipython3 # Baseline structural matrices for the true economy f = 1.05 @@ -405,7 +415,7 @@ The response shows the investment accelerator clearly: the full impact on net income $y_n$ occurs at lag 0, while consumption adjusts by only $1 - f^{-1} \approx 0.048$ and investment absorbs the remainder. -From lag 1 onward the economy is in its new steady state +From lag 1 onward the economy is in its new steady state, with permanently higher capital, income, and consumption ```{code-cell} ipython3 def table2_irf(A, C, n_lags=6): @@ -439,14 +449,14 @@ v_{t+1} = D v_t + \eta_t, where $\eta_t$ is a vector white noise with $E \eta_t \eta_t^\top = \Sigma_\eta$ and -$E \varepsilon_t v_s^\top = 0$ for all $t, s$. +$E w_t v_s^\top = 0$ for all $t, s$. The parameters are ```{math} -D = \operatorname{diag}(0.6, 0.7, 0.3), +D = \operatorname{diag}(\rho_y, \rho_c, \rho_{\Delta k}) = \operatorname{diag}(0.6, 0.7, 0.3), \qquad -\sigma_\eta = (0.05, 0.035, 0.65), +\Sigma_\eta = \operatorname{diag}(\sigma_{\eta}^2), \quad \sigma_\eta = (0.05, 0.035, 0.65), ``` so the unconditional covariance of $v_t$ is @@ -455,14 +465,15 @@ so the unconditional covariance of $v_t$ is R = \operatorname{diag}\!\left(\frac{\sigma_{\eta,i}^2}{1 - \rho_i^2}\right). ``` -The innovation variances are smallest for consumption -($\sigma_\eta = 0.035$), next for income ($\sigma_\eta = 0.05$), -and largest for investment ($\sigma_\eta = 0.65$). +The innovation standard deviations $\sigma_\eta$ are smallest for consumption +($0.035$), next for income ($0.05$), +and largest for investment ($0.65$). + +What matters for Granger-causality asymmetries, however, is the size of each measurement error **relative to the signal** in the corresponding true series. -As in {cite:t}`Sargent1989` and our discussion above, what matters for Granger-causality -asymmetries is the overall measurement quality in the full system: -output is relatively well measured while investment is relatively -poorly measured. +The innovations to the true series are $\theta_t$ for income (standard deviation $1$), $(1-\beta)\theta_t$ for consumption (standard deviation $0.048$), and $\beta \theta_t$ for investment (standard deviation $0.95$). + +Relative to these, income is very well measured (a ratio of about $0.05$), while consumption and investment are both poorly measured (ratios of about $0.7$). ```{code-cell} ipython3 ρ = np.array([0.6, 0.7, 0.3]) @@ -548,7 +559,7 @@ autoregressions for $\bar z_t$. Constructing the vector autoregression is also useful as an intermediate step in computing the likelihood of a sample of $\bar z_t$'s as a function of the free parameters -$\{A, C, D, Q, R\}$. +$\{A, C, D, Q, \Sigma_\eta\}$. The particular method that will be used to construct the vector autoregressive representation also proves useful as an intermediate @@ -576,7 +587,7 @@ Define ```{math} :label: model1_qd \tilde z_t = \bar z_{t+1} - D \bar z_t, \qquad -\bar\nu_t = C \varepsilon_t + \eta_t, \qquad +\bar\nu_t = C w_t + \eta_t, \qquad \bar C = CA - DC. ``` @@ -587,22 +598,24 @@ imply the state-space system ```{math} :label: model1_transformed \begin{aligned} -x_{t+1} &= A x_t + \varepsilon_t, \\ +x_{t+1} &= A x_t + w_t, \\ \tilde z_t &= \bar C\, x_t + \bar\nu_t, \end{aligned} ``` -where $(\varepsilon_t, \bar\nu_t)$ is a white noise process with +where $(w_t, \bar\nu_t)$ is a white noise process with ```{math} :label: model1_covs -E \begin{bmatrix} \varepsilon_t \end{bmatrix} -\begin{bmatrix} \varepsilon_t^\top & \bar\nu_t^\top \end{bmatrix} +E \begin{bmatrix} w_t \\ \bar\nu_t \end{bmatrix} +\begin{bmatrix} w_t \\ \bar\nu_t \end{bmatrix}^\top = \begin{bmatrix} Q & W_1 \\ W_1^\top & R_1 \end{bmatrix}, \qquad -R_1 = C Q C^\top + R, \quad W_1 = Q C^\top. +R_1 = C Q C^\top + \Sigma_\eta, \quad W_1 = Q C^\top. ``` +Notice that $R_1$ involves the covariance matrix $\Sigma_\eta$ of the **innovations** $\eta_t$ to the measurement errors, not the unconditional covariance matrix $R$ of the measurement errors themselves: quasi-differencing removes the serially correlated part of $v_t$. + System {eq}`model1_transformed` with covariances {eq}`model1_covs` is characterized by the five matrices $[A, \bar C, Q, R_1, W_1]$. @@ -667,7 +680,7 @@ Below we compute $K_1$, $S_1$, and $V_1$ numerically ```{code-cell} ipython3 C_bar = C @ A - D @ C -R1 = C @ Q @ C.T + R +R1 = C @ Q @ C.T + Σ_η W1 = Q @ C.T K1, S1, V1 = steady_state_kalman(A, C_bar, Q, R1, W1) @@ -686,9 +699,9 @@ r_t = \begin{bmatrix} \hat x_{t-1} \\ \bar z_{t-1} \end{bmatrix}, with dynamics ```{math} -r_{t+1} = F_1 r_t + G_1 u_t, +r_{t+1} = F_1 r_t + G_1 u_{t-1}, \qquad -\bar z_t = H_1 r_t + u_t, +\bar z_t = H_1 r_t + u_{t-1}, ``` where @@ -709,6 +722,8 @@ I H_1 = [\bar C \;\; D]. ``` +(Recall that $u_{t-1}$ is the innovation to $\bar z_t$.) + The moving average coefficients are then $\psi_0 = I$ and $\psi_j = H_1 F_1^{j-1} G_1$ for $j \geq 1$. @@ -759,12 +774,12 @@ $\hat x_0$, can be represented as ```{math} :label: model1_loglik -\mathcal{L}^* = -T\ln 2\pi - \tfrac{1}{2}T\ln|V_1| +\mathcal{L}^* = -\tfrac{3}{2}T\ln 2\pi - \tfrac{1}{2}T\ln|V_1| - \tfrac{1}{2}\sum_{t=0}^{T-1} u_t^\top V_1^{-1} u_t, ``` where $u_t$ is a function of $\{\bar z_t\}$ defined by -{eq}`model1_recursion` below. +{eq}`model1_recursion` below and the $3$ in the first term is the number of measured series. To use {eq}`model1_innov` to compute $\{u_t\}$, it is useful to represent it as @@ -792,7 +807,7 @@ $\{\bar z_t\}$. To measure the relative importance of each innovation, we decompose the $j$-step-ahead forecast-error variance of each measured variable. -Write $\bar z_{t+j} - E_t \bar z_{t+j} = \sum_{i=0}^{j-1} \psi_i u_{t+j-i}$. +Write $\bar z_{t+j} - E_t \bar z_{t+j} = \sum_{i=0}^{j-1} \psi_i u_{t+j-i}$, where, to lighten notation, here $u_{t+j-i}$ denotes the innovation to $\bar z_{t+j-i}$. Let $P$ be the lower-triangular Cholesky factor of $V_1$ so that the orthogonalized innovations are $e_t = P^{-1} u_t$. @@ -801,6 +816,10 @@ Then the contribution of orthogonalized innovation $k$ to the $j$-step-ahead variance of variable $m$ is $\sum_{i=0}^{j-1} (\psi_i P)_{mk}^2$. +We order the variables as $y_n$, $c$, $\Delta k$. + +With this ordering, the first orthogonalized innovation is the innovation in $y_n$, the second is the part of the innovation in $c$ that is orthogonal to the innovation in $y_n$, and the third is the part of the innovation in $\Delta k$ that is orthogonal to the first two. + The table below shows the cumulative contribution of each orthogonalized innovation to the forecast-error variance of $y_n$, $c$, and $\Delta k$ at horizons 1 through 20. @@ -846,6 +865,10 @@ Granger-cause consumption and investment, but not vice versa. This matches the paper's message that, in a one-common-index model, the relatively best measured series has the strongest predictive content. +Because income is ordered first, the decomposition attributes to the income innovation all of the variance that it shares with the other innovations. + +Below we confirm the Granger-causality pattern with a calculation that does not depend on an ordering. + Let's look at the covariance matrix of the innovations ```{code-cell} ipython3 @@ -866,6 +889,54 @@ print(np.sort(np.linalg.eigvalsh(V1))[::-1].round(4)) The first eigenvalue is much larger than the others, consistent with the presence of a dominant common shock $\theta_t$ +### Granger causality without an ordering + +The population vector autoregression for $\bar z_t$ gives a direct view of Granger causality. + +Iterating on {eq}`model1_recursion` gives + +$$ +\tilde z_t = \sum_{j=1}^\infty \bar C (A - K_1 \bar C)^{j-1} K_1 \tilde z_{t-j} + u_t . +$$ + +Let $B_j = \bar C (A - K_1 \bar C)^{j-1} K_1$. + +Substituting $\tilde z_t = \bar z_{t+1} - D \bar z_t$ and collecting terms gives the vector autoregression + +$$ +\bar z_{t+1} = \sum_{j=1}^\infty \Pi_j \bar z_{t+1-j} + u_t, +\qquad +\Pi_1 = D + B_1, \quad \Pi_j = B_j - B_{j-1} D \ \ (j \geq 2). +$$ + +Variable $i$ fails to Granger-cause variable $m$ if and only if the $(m, i)$ entries of all the $\Pi_j$ are zero. + +Here are the first two coefficient matrices; rows are equations and columns are lagged variables. + +```{code-cell} ipython3 +M_filter = A - K1 @ C_bar +print("eigenvalues of A - K1 C_bar:", np.linalg.eigvals(M_filter).round(4)) + +B_coef = [C_bar @ np.linalg.matrix_power(M_filter, j) @ K1 for j in range(3)] +Π = [D + B_coef[0]] + [B_coef[j] - B_coef[j-1] @ D for j in range(1, 3)] + +eq_labels = [r'\text{eq. } ' + lab for lab in labels] +for j in range(2): + df = pd.DataFrame(Π[j], index=eq_labels, columns=labels) + display(Latex(r'$\Pi_' + str(j + 1) + ' = ' + + df_to_latex_array(df).strip('$') + '$')) +``` + +The eigenvalues of $A - K_1 \bar C$ are inside the unit circle, so the vector autoregression converges. + +In the investment equation, the coefficient on lagged income is large (about $-0.28$). + +In the income equation, the coefficient on lagged investment is tiny (about $-0.004$). + +Lagged income helps to predict investment, but lagged investment barely helps to predict income, which is the Granger-causality pattern that {cite:t}`Sargent1989` set out to explain. + +{ref}`meas_ex2` computes an ordering-free measure of the strength of Granger causality due to {cite:t}`Geweke1982`. + ### Wold impulse responses Impulse responses in the Wold representation are reported using orthogonalized @@ -913,7 +984,7 @@ series, its innovation is dominated by the true permanent shock $\theta_t$. The consumption and investment innovations produce responses that -decay according to the AR(1) structure of their respective +decay at rates close to the autoregressive coefficients of their respective measurement errors ($\rho_c = 0.7$, $\rho_{\Delta k} = 0.3$), with little spillover to other variables. @@ -948,10 +1019,10 @@ To prepare its estimates, the reporting agency itself computes the Kalman filter to obtain the innovations representation {eq}`model1_innov`. Rather than reporting the error-corrupted data $\bar z_t$, the agency -reports $\tilde z_t = G \hat x_t$, where $G$ is a "selection matrix," +reports $\hat z_t = G \hat x_t$, where $G$ is a "selection matrix," possibly equal to $C$, for the data reported by the agency. -The data $G \hat x_t = E[G x_t \mid \bar z_t, \bar z_{t-1}, \ldots, \hat x_0]$. +The reported data are $\hat z_t = G \hat x_t = E[G x_t \mid \bar z_t, \bar z_{t-1}, \ldots, \hat x_0]$. The state-space representation for the reported data is then @@ -959,7 +1030,7 @@ The state-space representation for the reported data is then :label: model2_state \begin{aligned} \hat x_{t+1} &= A \hat x_t + K_1 u_t, \\ -\tilde z_t &= G \hat x_t, +\hat z_t &= G \hat x_t, \end{aligned} ``` @@ -967,10 +1038,10 @@ where the first line of {eq}`model2_state` is from the innovations representation {eq}`model1_innov`. Note that $u_t$ is the innovation to $\bar z_{t+1}$ and is *not* the -innovation to $\tilde z_t$. +innovation to $\hat z_t$. -To obtain a Wold representation for $\tilde z_t$ and the likelihood -function for a sample of $\tilde z_t$ requires that we obtain an +To obtain a Wold representation for $\hat z_t$ and the likelihood +function for a sample of $\hat z_t$ requires that we obtain an innovations representation for {eq}`model2_state`. ### Innovations representation for filtered data @@ -981,19 +1052,19 @@ To add a little generality to {eq}`model2_state` we amend it to the system :label: model2_obs \begin{aligned} \hat x_{t+1} &= A \hat x_t + K_1 u_t, \\ -\tilde z_t &= G \hat x_t + \eta_t, +\hat z_t &= G \hat x_t + \zeta_t, \end{aligned} ``` -where $\eta_t$ is a type 2 white-noise measurement error process +where $\zeta_t$ is a type 2 white-noise measurement error process ("typos") with presumably very small covariance matrix $R_2$. The covariance matrix of the joint noise is ```{math} :label: model2_Q -E \begin{bmatrix} K_1 u_t \\ \eta_t \end{bmatrix} - \begin{bmatrix} K_1 u_t \\ \eta_t \end{bmatrix}^\top +E \begin{bmatrix} K_1 u_t \\ \zeta_t \end{bmatrix} + \begin{bmatrix} K_1 u_t \\ \zeta_t \end{bmatrix}^\top = \begin{bmatrix} Q_2 & 0 \\ 0 & R_2 \end{bmatrix}, ``` @@ -1013,7 +1084,7 @@ representation is :label: model2_innov \begin{aligned} \check{x}_{t+1} &= A \check{x}_t + K_2 a_t, \\ -\tilde z_t &= G \check{x}_t + a_t, +\hat z_t &= G \check{x}_t + a_t, \end{aligned} ``` @@ -1022,15 +1093,15 @@ where ```{math} :label: model2_innov_defs \begin{aligned} -a_t &= \tilde z_t - E[\tilde z_t \mid \tilde z_{t-1}, \tilde z_{t-2}, \ldots], \\ -\check{x}_t &= E[\hat x_t \mid \tilde z_{t-1}, \tilde z_{t-2}, \ldots, \check{x}_0], \\ +a_t &= \hat z_t - E[\hat z_t \mid \hat z_{t-1}, \hat z_{t-2}, \ldots], \\ +\check{x}_t &= E[\hat x_t \mid \hat z_{t-1}, \hat z_{t-2}, \ldots, \check{x}_0], \\ S_2 &= E[(\hat x_t - \check{x}_t)(\hat x_t - \check{x}_t)^\top], \\ [K_2, S_2] &= \text{kalmanfilter}(A, G, Q_2, R_2, 0). \end{aligned} ``` Thus $\{a_t\}$ is the innovation process for the reported data -$\tilde z_t$, with innovation covariance +$\hat z_t$, with innovation covariance ```{math} :label: model2_V2 @@ -1039,12 +1110,12 @@ V_2 = E\, a_t a_t^\top = G\, S_2\, G^\top + R_2. ### Wold representation -A Wold moving average representation for $\tilde z_t$ is found from +A Wold moving average representation for $\hat z_t$ is found from {eq}`model2_innov` to be ```{math} :label: model2_wold -\tilde z_t = \bigl[G(I - AL)^{-1} K_2 L + I\bigr] a_t, +\hat z_t = \bigl[G(I - AL)^{-1} K_2 L + I\bigr] a_t, ``` with coefficients $\psi_0 = I$ and $\psi_j = G A^{j-1} K_2$ for @@ -1057,23 +1128,23 @@ to undo. ### Gaussian likelihood When a method analogous to Model 1 is used, a Gaussian log-likelihood -for $\tilde z_t$ can be computed by first computing an $\{a_t\}$ sequence -from observations on $\tilde z_t$ by using +for $\hat z_t$ can be computed by first computing an $\{a_t\}$ sequence +from observations on $\hat z_t$ by using ```{math} :label: model2_recursion \begin{aligned} -\check{x}_{t+1} &= (A - K_2 G)\,\check{x}_t + K_2 \tilde z_t, \\ -a_t &= -G\,\check{x}_t + \tilde z_t. +\check{x}_{t+1} &= (A - K_2 G)\,\check{x}_t + K_2 \hat z_t, \\ +a_t &= -G\,\check{x}_t + \hat z_t. \end{aligned} ``` The likelihood function for a sample of $T$ observations -$\{\tilde z_t\}$ is then +$\{\hat z_t\}$ is then ```{math} :label: model2_loglik -\mathcal{L}^{**} = -T\ln 2\pi - \tfrac{1}{2}T\ln|V_2| +\mathcal{L}^{**} = -\tfrac{3}{2}T\ln 2\pi - \tfrac{1}{2}T\ln|V_2| - \tfrac{1}{2}\sum_{t=0}^{T-1} a_t^\top V_2^{-1} a_t. ``` @@ -1099,6 +1170,10 @@ it does not need the Wold representation for the filtered data. In our parameterization $G = C$. +Because the rows of $C$ satisfy the national income identity, the reported data satisfy $\hat c_t + \Delta \hat k_t = \hat y_{nt}$ exactly. + +So when $R_2 = 0$, the innovation covariance matrix $V_2$ is singular, and the likelihood {eq}`model2_loglik` of all three reported series is degenerate; {ref}`meas_ex3` explores this "stochastic singularity." + ```{code-cell} ipython3 Q2 = K1 @ V1 @ K1.T ε = 1e-6 @@ -1161,6 +1236,12 @@ print('Eigenvalues of covariance matrix:') print(np.sort(np.linalg.eigvalsh(V2))[::-1].round(4)) ``` +One eigenvalue is large. + +A second eigenvalue is positive but tiny: it reflects revisions to the agency's estimate of the capital stock, which enter the reported data only through the small coefficient $f - 1$. + +The third eigenvalue equals $\epsilon$ and would be exactly zero if $R_2 = 0$, because of the national income identity. + As {cite:t}`Sargent1989` emphasizes, the two models of measurement produce quite different inferences about the economy's dynamics despite sharing identical underlying parameters. @@ -1250,7 +1331,7 @@ def simulate_series(seed=7909, T=80, k0=10.0): # Filtered data via Model 1 transformed filter xhat_prev = np.array([k0, 0.0]) - z_prev = np.zeros(3) + z_prev = C @ xhat_prev # model-implied value of the unobserved z̄_{-1} z_filt = np.empty((T, 3)) k_filt = np.empty(T) @@ -1398,13 +1479,15 @@ plot_true_vs_other(t, sim["k_true"], sim["k_filt"], Indeed, Kalman-filtered estimates from Model 1 remove much of the measurement noise and track the truth closely. +The improvement is largest for investment, the series whose measurement error is largest; income was already well measured. + In the true model the national income identity $c_t + \Delta k_t = y_{n,t}$ holds exactly. Independent measurement errors break this accounting identity in the measured data. -The Kalman filter approximately restores it. +Because the filtered series are all of the form $C \hat x_t$, and the rows of $C$ satisfy the identity, the filtered data satisfy it exactly. The following figure confirms this by showing the residual $c_t + \Delta k_t - y_{n,t}$ for both measured and filtered data @@ -1432,7 +1515,7 @@ plt.tight_layout() plt.show() ``` -As we have predicted, the residual for the measured data is large and volatile, while the residual for the filtered data is numerically 0. +As we have predicted, the residual for the measured data is large and volatile, while the residual for the filtered data is zero up to rounding error. ## Summary @@ -1452,6 +1535,8 @@ the empirical accelerator literature: income appears to Granger-cause consumption and investment, a result {cite:t}`Sargent1989` attributes to measurement error and signal extraction in raw reported data. +What matters for this pattern is how large each series' measurement error is relative to the variation in the true series. + Model 2, working with filtered data, attributes nearly all variance to the single structural shock $\theta_t$ and *cannot* reproduce the Granger causality pattern. @@ -1460,5 +1545,179 @@ The {doc}`Kalman filter ` effectively strips measurement noise from the data, so the filtered series track the truth closely. Raw measurement error breaks the national income accounting identity, -but the near-zero residual shows that the filter approximately -restores it. +while the filtered data satisfy it exactly, which makes the three reported series stochastically singular. + +## Exercises + +```{exercise} +:label: meas_ex1 + +Under the assumptions $a(L) = 1$ and $\beta f = 1$, use the decision rule {eq}`simple_crule` and the law of motion $k_{t+1} = k_t + f^{-1} \theta_t$ to show that + +1. $c_t - c_{t-1} = (1 - \beta)\theta_t$, so that consumption is a martingale +1. $c_t - \beta c_{t-1} = (1-\beta) y_{nt}$, which is Friedman's consumption function {eq}`friedman_consumption` +1. $y_{nt} - y_{n,t-1} = \theta_t - \beta \theta_{t-1}$, which is equation {eq}`income_process` +1. $(k_{t+1} - k_t) - \beta (k_t - k_{t-1}) = \beta (y_{nt} - y_{n,t-1})$, which is the accelerator {eq}`mm_accelerator` + +Then verify all four identities numerically using the true series in `sim`, recovering $\theta_t$ from $\theta_t = y_{nt} - (f-1)k_t$. +``` + +```{solution-start} meas_ex1 +:class: dropdown +``` + +From $k_t - k_{t-1} = \beta \theta_{t-1}$ and {eq}`simple_crule`, + +$$ +c_t - c_{t-1} = (f-1)\beta\theta_{t-1} + (1-\beta)(\theta_t - \theta_{t-1}) = (1-\beta)\theta_t +$$ + +because $(f-1)\beta = 1 - \beta$. + +Similarly, $y_{nt} - y_{n,t-1} = (f-1)\beta \theta_{t-1} + \theta_t - \theta_{t-1} = \theta_t - \beta\theta_{t-1}$. + +For the third claim, substitute $k_t = k_{t-1} + \beta\theta_{t-1}$ into {eq}`simple_crule` and {eq}`net_income`: + +$$ +\begin{aligned} +c_t - \beta c_{t-1} & = (f-1)(1-\beta) k_{t-1} + (1-\beta)\theta_t + \left[(f-1)\beta - \beta(1-\beta)\right]\theta_{t-1} \\ +(1-\beta) y_{nt} & = (1-\beta)(f-1) k_{t-1} + (1-\beta)\theta_t + (1-\beta)(f-1)\beta\,\theta_{t-1} +\end{aligned} +$$ + +and use $(f-1)\beta = 1 - \beta$ to see that the two right sides are equal. + +Finally, $k_{t+1} - k_t = \beta\theta_t$ implies $(1 - \beta L)(k_{t+1} - k_t) = \beta(\theta_t - \beta\theta_{t-1}) = \beta (y_{nt} - y_{n,t-1})$. + +```{code-cell} ipython3 +y, c, dk, k = sim["y_true"], sim["c_true"], sim["dk_true"], sim["k_true"] +θ = y - (f - 1) * k + +print("martingale: ", np.allclose(np.diff(c), (1 - β) * θ[1:])) +print("Friedman: ", np.allclose(c[1:] - β * c[:-1], (1 - β) * y[1:])) +print("MA(1): ", np.allclose(np.diff(y), θ[1:] - β * θ[:-1])) +print("accelerator: ", np.allclose(dk[1:] - β * dk[:-1], β * np.diff(y))) +``` + +```{solution-end} +``` + +```{exercise} +:label: meas_ex2 + +Forecast-error-variance decompositions based on a Cholesky factorization depend on how the variables are ordered. + +{cite:t}`Geweke1982` proposed a measure of Granger causality that does not depend on an ordering. + +For measured series $i$ and $j$, let $\sigma^2_{j}$ be the one-step-ahead prediction error variance of $\bar z_{j,t+1}$ given its own past only, and let $\sigma^2_{j \mid i}$ be its prediction error variance given the past of both $\bar z_{i}$ and $\bar z_{j}$. + +Then $F_{i \to j} = \log(\sigma^2_j / \sigma^2_{j \mid i}) \geq 0$, with equality if and only if $\bar z_i$ does not Granger-cause $\bar z_j$. + +Because the measurement errors are mutually independent AR(1) processes, both variances can be computed by applying `steady_state_kalman` to the quasi-differenced system {eq}`model1_transformed` using only the rows of $\bar C$, $R_1$, and the columns of $W_1$ that correspond to the included series. + +1. Compute $F_{y \to \Delta k}$, $F_{\Delta k \to y}$, $F_{y \to c}$, and $F_{c \to y}$ for the baseline measurement-error parameters. +1. Recompute them after swapping the measurement-error innovation standard deviations of income and investment, so that $\sigma_\eta = (0.65, 0.035, 0.05)$. +1. Recompute them with $\sigma_\eta = (0.3, 0.3, 0.3)$ and explain why equal measurement error variances do **not** produce symmetric Granger causality. +``` + +```{solution-start} meas_ex2 +:class: dropdown +``` + +Here is one solution. + +```{code-cell} ipython3 +def geweke(σ_η, ρ=np.array([0.6, 0.7, 0.3]), pair=(0, 2)): + """ + Return (F_{i->j}, F_{j->i}) for the measured series pair = (i, j). + """ + D = np.diag(ρ) + Σ = np.diag(np.asarray(σ_η)**2) + C_bar = C @ A - D @ C + R1 = C @ Q @ C.T + Σ + W1 = Q @ C.T + + def pev(idx): + idx = list(idx) + return steady_state_kalman(A, C_bar[idx], Q, + R1[np.ix_(idx, idx)], W1[:, idx])[2] + + i, j = pair + V_i, V_j, V_ij = pev([i])[0, 0], pev([j])[0, 0], pev([i, j]) + return np.log(V_j / V_ij[1, 1]), np.log(V_i / V_ij[0, 0]) + + +cases = {"baseline": (0.05, 0.035, 0.65), + "swapped y and Δk": (0.65, 0.035, 0.05), + "equal": (0.3, 0.3, 0.3)} + +for name, σ in cases.items(): + F_yk, F_ky = geweke(σ, pair=(0, 2)) + F_yc, F_cy = geweke(σ, pair=(0, 1)) + print(f"{name:18s} y→Δk {F_yk:.4f} Δk→y {F_ky:.4f} " + f"y→c {F_yc:.4f} c→y {F_cy:.4f}") +``` + +In the baseline, income Granger-causes both investment and consumption, while neither Granger-causes income. + +Swapping the measurement errors of income and investment reverses the direction of Granger causality between them. + +With equal measurement error variances, investment Granger-causes income more than income Granger-causes investment, and income strongly Granger-causes consumption. + +What matters is not the variance of the measurement error but its size **relative to the signal** in each series. + +The innovation to true investment $\beta \theta_t$ has standard deviation $0.95$, while the innovation to true consumption $(1-\beta)\theta_t$ has standard deviation only $0.048$. + +So a measurement error with innovation standard deviation $0.3$ obscures investment much less than it swamps consumption. + +In the baseline, consumption's measurement error innovation ($0.035$) is small in absolute terms but not relative to the tiny innovation in true consumption, which is why income Granger-causes consumption. + +```{solution-end} +``` + +```{exercise} +:label: meas_ex3 + +In Model 2, the agency reports $\hat z_t = C \hat x_t$. + +1. Show that the reported data satisfy the national income identity $\hat c_t + \Delta \hat k_t - \hat y_{nt} = 0$ **exactly**, for any estimate $\hat x_t$. +1. Explain why this implies that the innovation covariance matrix $V_2$ of the reported data is singular when $R_2 = 0$. +1. Compute the eigenvalues of $V_2$ for $R_2 = \epsilon I$ with $\epsilon = 10^{-2}, 10^{-4}, 10^{-6}$ and interpret them. +1. What does this "stochastic singularity" imply for an econometrician who wants to compute the likelihood {eq}`model2_loglik` from all three reported series? +``` + +```{solution-start} meas_ex3 +:class: dropdown +``` + +The rows of $C$ satisfy $C_{c} + C_{\Delta k} - C_{y} = [(f-1) + 0 - (f-1), \; (1 - f^{-1}) + f^{-1} - 1] = [0, 0]$. + +So $\hat c_t + \Delta \hat k_t - \hat y_{nt} = (C_c + C_{\Delta k} - C_y)\hat x_t = 0$ for every $\hat x_t$. + +Because the vector $b = [-1, 1, 1]$ satisfies $b\, C \hat x_t = 0$ for all $t$, it also satisfies $b\, a_t = 0$ when $R_2 = 0$, so $b V_2 b^\top = 0$ and $V_2$ is singular. + +```{code-cell} ipython3 +print("C_c + C_Δk - C_y =", C[1] + C[2] - C[0]) +print("max |identity residual| in reported data:", + np.max(np.abs(sim['c_filt'] + sim['dk_filt'] - sim['y_filt']))) + +for ε_R in (1e-2, 1e-4, 1e-6): + _, _, V2_ε = steady_state_kalman(A, C, Q2, ε_R * np.eye(3)) + print(f"ε = {ε_R:.0e}: eigenvalues of V2 = " + f"{np.sort(np.linalg.eigvalsh(V2_ε))[::-1]}") +``` + +One eigenvalue of $V_2$ is large: it reflects news about the permanent shock $\theta_t$. + +A second eigenvalue is tiny but positive: it reflects revisions to the agency's estimate of the capital stock, which enter the reported data only through the small coefficient $f - 1$. + +The third eigenvalue is approximately $\epsilon$ and vanishes as $\epsilon \to 0$ because of the identity. + +An econometrician who tries to compute the Gaussian likelihood of all three reported series with $R_2 = 0$ faces a degenerate distribution: $\ln |V_2| = -\infty$. + +In practice, one must either drop one of the three series, or add "type 2" reporting errors with $R_2 > 0$ as in {eq}`model2_obs`. + +In the raw data of Model 1, by contrast, independent measurement errors break the identity, so the three measured series have a nonsingular innovation covariance matrix. + +```{solution-end} +``` diff --git a/lectures/svd_intro.md b/lectures/svd_intro.md index daae7edc6..4aca757af 100644 --- a/lectures/svd_intro.md +++ b/lectures/svd_intro.md @@ -16,47 +16,48 @@ kernelspec: ## Overview The **singular value decomposition** (SVD) is a work-horse in applications of least squares projection that -form foundations for many statistical and machine learning methods. +form foundations for many statistical and machine learning methods. After defining the SVD, we'll describe how it connects to * **four fundamental spaces** of linear algebra -* under-determined and over-determined **least squares regressions** +* **low-rank approximations** of matrices +* under-determined and over-determined **least squares regressions** (see the exercises) * **principal components analysis** (PCA) -Like principal components analysis (PCA), DMD can be thought of as a data-reduction procedure that represents salient patterns by projecting data onto a limited set of factors. +In a sequel to this lecture, {doc}`VARs and DMDs `, we'll describe how SVDs provide ways rapidly to compute reduced-rank approximations to first-order vector autoregressions (VARs) by means of a **dynamic mode decomposition** (DMD). -In a sequel to this lecture about {doc}`Dynamic Mode Decompositions `, we'll describe how SVD's provide ways rapidly to compute reduced-order approximations to first-order Vector Autoregressions (VARs). +Like principal components analysis, DMD can be thought of as a data-reduction procedure that represents salient patterns by projecting data onto a limited set of factors. -## The Setting +## The Setting Let $X$ be an $m \times n$ matrix of rank $p$. Necessarily, $p \leq \min(m,n)$. -In much of this lecture, we'll think of $X$ as a matrix of data in which +In much of this lecture, we'll think of $X$ as a matrix of data in which * each column is an **individual** -- a time period or person, depending on the application * each row is a **random variable** describing an attribute of a time period or a person, depending on the application -We'll be interested in two situations +We'll be interested in two situations -* A **short and fat** case in which $m << n$, so that there are many more columns (individuals) than rows (attributes). +* A **short and fat** case in which $m \ll n$, so that there are many more columns (individuals) than rows (attributes). -* A **tall and skinny** case in which $m >> n$, so that there are many more rows (attributes) than columns (individuals). +* A **tall and skinny** case in which $m \gg n$, so that there are many more rows (attributes) than columns (individuals). We'll apply a **singular value decomposition** of $X$ in both situations. -In the $ m < < n$ case in which there are many more individuals $n$ than attributes $m$, we can calculate sample moments of a joint distribution by taking averages across observations of functions of the observations. +In the $m \ll n$ case in which there are many more individuals $n$ than attributes $m$, we can calculate sample moments of a joint distribution by taking averages across observations of functions of the observations. -In this $ m < < n$ case, we'll look for patterns by using a singular value decomposition to do a principal components analysis (PCA). +In this $m \ll n$ case, we'll look for patterns by using a singular value decomposition to do a principal components analysis (PCA). -In the $m > > n$ case in which there are many more attributes $m$ than individuals $n$ and when we are in a time-series setting in which $n$ equals the number of time periods covered in the data set $X$, we'll proceed in a different way. +In the $m \gg n$ case in which there are many more attributes $m$ than individuals $n$ and when we are in a time-series setting in which $n$ equals the number of time periods covered in the data set $X$, we'll proceed in a different way. -We'll again use a singular value decomposition, but now to construct a **dynamic mode decomposition** (DMD) +In the sequel {doc}`VARs and DMDs `, we'll again use a singular value decomposition, but now to construct a **dynamic mode decomposition** (DMD). ## Singular Value Decomposition @@ -77,13 +78,13 @@ $$ and -* $U$ is an $m \times m$ orthogonal matrix of **left singular vectors** of $X$ +* $U$ is an $m \times m$ orthogonal matrix of **left singular vectors** of $X$ * Columns of $U$ are eigenvectors of $X X^\top $ * $V$ is an $n \times n$ orthogonal matrix of **right singular vectors** of $X$ -* Columns of $V$ are eigenvectors of $X^\top X$ -* $\Sigma$ is an $m \times n$ matrix in which the first $p$ places on its main diagonal are positive numbers $\sigma_1, \sigma_2, \ldots, \sigma_p$ called **singular values**; remaining entries of $\Sigma$ are all zero +* Columns of $V$ are eigenvectors of $X^\top X$ +* $\Sigma$ is an $m \times n$ matrix in which the first $p$ places on its main diagonal are positive numbers $\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_p$ called **singular values**; remaining entries of $\Sigma$ are all zero -* The $p$ singular values are positive square roots of the eigenvalues of the $m \times m$ matrix $X X^\top $ and also of the $n \times n$ matrix $X^\top X$ +* The $p$ singular values are positive square roots of the nonzero eigenvalues of the $m \times m$ matrix $X X^\top $ and also of the $n \times n$ matrix $X^\top X$ * We adopt a convention that when $U$ is a complex valued matrix, $U^\top $ denotes the **conjugate-transpose** or **Hermitian-transpose** of $U$, meaning that $U_{ij}^\top $ is the complex conjugate of $U_{ji}$. @@ -91,34 +92,34 @@ $U_{ij}^\top $ is the complex conjugate of $U_{ji}$. * Similarly, when $V$ is a complex valued matrix, $V^\top $ denotes the **conjugate-transpose** or **Hermitian-transpose** of $V$ -The matrices $U,\Sigma,V$ entail linear transformations that reshape in vectors in the following ways: +The matrices $U,\Sigma,V$ entail linear transformations that reshape vectors in the following ways: -* multiplying vectors by the unitary matrices $U$ and $V$ **rotates** them, but leaves **angles between vectors** and **lengths of vectors** unchanged. -* multiplying vectors by the diagonal matrix $\Sigma$ leaves **angles between vectors** unchanged but **rescales** vectors. +* multiplying vectors by the unitary matrices $U$ and $V$ **rotates** them, but leaves **angles between vectors** and **lengths of vectors** unchanged. +* multiplying vectors by the diagonal matrix $\Sigma$ **rescales** each coordinate by a singular value (and adds or drops coordinates when $m \neq n$); in general this changes both lengths of vectors and angles between them. -Thus, representation {eq}`eq:SVD101` asserts that multiplying an $n \times 1$ vector $y$ by the $m \times n$ matrix $X$ +Thus, representation {eq}`eq:SVD101` asserts that multiplying an $n \times 1$ vector $y$ by the $m \times n$ matrix $X$ amounts to performing the following three multiplications of $y$ sequentially: * **rotating** $y$ by computing $V^\top y$ * **rescaling** $V^\top y$ by multiplying it by $\Sigma$ * **rotating** $\Sigma V^\top y$ by multiplying it by $U$ -This structure of the $m \times n$ matrix $X$ opens the door to constructing systems +This structure of the $m \times n$ matrix $X$ opens the door to constructing systems of data **encoders** and **decoders**. Thus, -* $V^\top y$ is an encoder +* $V^\top$ is an encoder * $\Sigma$ is an operator to be applied to the encoded data * $U$ is a decoder to be applied to the output from applying operator $\Sigma$ to the encoded data -We'll apply this circle of ideas later in this lecture when we study Dynamic Mode Decomposition. +We'll apply this circle of ideas in the sequel {doc}`VARs and DMDs ` when we study dynamic mode decompositions. **Road Ahead** -What we have described above is called a **full** SVD. +What we have described above is called a **full** SVD. -In a **full** SVD, the shapes of $U$, $\Sigma$, and $V$ are $\left(m, m\right)$, $\left(m, n\right)$, $\left(n, n\right)$, respectively. +In a **full** SVD, the shapes of $U$, $\Sigma$, and $V$ are $\left(m, m\right)$, $\left(m, n\right)$, $\left(n, n\right)$, respectively. Later we'll also describe an **economy** or **reduced** SVD. @@ -126,20 +127,20 @@ Before we study a **reduced** SVD we'll say a little more about properties of a ## Four Fundamental Subspaces -Let ${\mathcal C}$ denote a column space, ${\mathcal N}$ denote a null space, and ${\mathcal R}$ denote a row space. +Let ${\mathcal C}$ denote a column space, ${\mathcal N}$ denote a null space, and ${\mathcal R}$ denote a row space. Let's start by recalling the four fundamental subspaces of an $m \times n$ matrix $X$ of rank $p$. -* The **column space** of $X$, denoted ${\mathcal C}(X)$, is the span of the columns of $X$, i.e., all vectors $y$ that can be written as linear combinations of columns of $X$. Its dimension is $p$. +* The **column space** of $X$, denoted ${\mathcal C}(X)$, is the span of the columns of $X$, i.e., all vectors $y$ that can be written as linear combinations of columns of $X$. Its dimension is $p$. * The **null space** of $X$, denoted ${\mathcal N}(X)$ consists of all vectors $y$ that satisfy $X y = 0$. Its dimension is $n-p$. * The **row space** of $X$, denoted ${\mathcal R}(X)$ is the column space of $X^\top $. It consists of all -vectors $z$ that can be written as linear combinations of rows of $X$. Its dimension is $p$. -* The **left null space** of $X$, denoted ${\mathcal N}(X^\top )$, consist of all vectors $z$ such that -$X^\top z =0$. Its dimension is $m-p$. +vectors $z$ that can be written as linear combinations of rows of $X$. Its dimension is $p$. +* The **left null space** of $X$, denoted ${\mathcal N}(X^\top )$, consists of all vectors $z$ such that +$X^\top z =0$. Its dimension is $m-p$. -For a full SVD of a matrix $X$, the matrix $U$ of left singular vectors and the matrix $V$ of right singular vectors contain orthogonal bases for all four subspaces. +For a full SVD of a matrix $X$, the matrix $U$ of left singular vectors and the matrix $V$ of right singular vectors contain orthogonal bases for all four subspaces. They form two pairs of orthogonal subspaces that we'll describe now. @@ -147,19 +148,19 @@ that we'll describe now. Let $u_i, i = 1, \ldots, m$ be the $m$ column vectors of $U$ and let $v_i, i = 1, \ldots, n$ be the $n$ column vectors of $V$. -Let's write the full SVD of X as +Let's write the full SVD of $X$ as $$ X = \begin{bmatrix} U_L & U_R \end{bmatrix} \begin{bmatrix} \Sigma_p & 0 \cr 0 & 0 \end{bmatrix} \begin{bmatrix} V_L & V_R \end{bmatrix}^\top $$ (eq:fullSVDpartition) -where $ \Sigma_p$ is a $p \times p$ diagonal matrix with the $p$ singular values on the diagonal and +where $\Sigma_p$ is a $p \times p$ diagonal matrix with the $p$ singular values on the diagonal and $$ \begin{aligned} -U_L & = \begin{bmatrix}u_1 & \cdots & u_p \end{bmatrix}, \quad U_R = \begin{bmatrix}u_{p+1} & \cdots u_m \end{bmatrix} \cr -V_L & = \begin{bmatrix}v_1 & \cdots & v_p \end{bmatrix} , \quad U_R = \begin{bmatrix}v_{p+1} & \cdots u_n \end{bmatrix} +U_L & = \begin{bmatrix}u_1 & \cdots & u_p \end{bmatrix}, \quad U_R = \begin{bmatrix}u_{p+1} & \cdots & u_m \end{bmatrix} \cr +V_L & = \begin{bmatrix}v_1 & \cdots & v_p \end{bmatrix} , \quad V_R = \begin{bmatrix}v_{p+1} & \cdots & v_n \end{bmatrix} \end{aligned} $$ @@ -188,7 +189,7 @@ X v_i & = 0 , \quad i = p+1, \ldots, n \end{aligned} $$ (eq:orthoortho1) -Equations {eq}`eq:orthoortho1` tell how the transformation $X$ maps a pair of orthonormal vectors $v_i, v_j$ for $i$ and $j$ both less than or equal to the rank $p$ of $X$ into a pair of orthonormal vectors $u_i, u_j$. +Equations {eq}`eq:orthoortho1` tell how the transformation $X$ maps a pair of orthonormal vectors $v_i, v_j$ for $i$ and $j$ both less than or equal to the rank $p$ of $X$ into a pair of orthogonal vectors $\sigma_i u_i, \sigma_j u_j$ with lengths $\sigma_i$ and $\sigma_j$. Equations {eq}`eq:Xfour1a` assert that @@ -225,7 +226,7 @@ X^\top u_i & = 0 \quad i= p+1, \ldots, m \end{aligned} $$ (eq:orthoortho2) -Notice how equations {eq}`eq:orthoortho2` assert that the transformation $X^\top $ maps a pair of distinct orthonormal vectors $u_i, u_j$ for $i$ and $j$ both less than or equal to the rank $p$ of $X$ into a pair of distinct orthonormal vectors $v_i, v_j$ . +Notice how equations {eq}`eq:orthoortho2` assert that the transformation $X^\top $ maps a pair of distinct orthonormal vectors $u_i, u_j$ for $i$ and $j$ both less than or equal to the rank $p$ of $X$ into a pair of orthogonal vectors $\sigma_i v_i, \sigma_j v_j$. Equations {eq}`eq:Xfour1b` assert that @@ -240,15 +241,14 @@ $$ Thus, taken together, the systems of equations {eq}`eq:Xfour1a` and {eq}`eq:Xfour1b` -describe the four fundamental subspaces of $X$ in the following ways: +describe the four fundamental subspaces of $X$ in the following ways: $$ \begin{aligned} {\mathcal C}(X) & = {\mathcal C}(U_L) \cr {\mathcal N}(X^\top ) & = {\mathcal C}(U_R) \cr {\mathcal R}(X) & \equiv {\mathcal C}(X^\top ) = {\mathcal C} (V_L) \cr -{\mathcal N}(X) & = {\mathcal C} (V_R) \cr - +{\mathcal N}(X) & = {\mathcal C} (V_R) \end{aligned} $$ (eq:fourspaceSVD) @@ -260,9 +260,9 @@ Since $U$ and $V$ are both orthonormal matrices, collection {eq}`eq:fourspaceSVD * $V_R$ is an orthonormal basis for the null space of $X$ -We have verified the four claims in {eq}`eq:fourspaceSVD` simply by performing the multiplications called for by the right side of {eq}`eq:fullSVDpartition` and reading them. +We have verified the four claims in {eq}`eq:fourspaceSVD` simply by performing the multiplications called for by the right side of {eq}`eq:fullSVDpartition` and reading them. -The claims in {eq}`eq:fourspaceSVD` and the fact that $U$ and $V$ are both unitary (i.e, orthonormal) matrices imply +The claims in {eq}`eq:fourspaceSVD` and the fact that $U$ and $V$ are both unitary (i.e., orthonormal) matrices imply that * the column space of $X$ is orthogonal to the null space of $X^\top $ @@ -270,8 +270,8 @@ that Sometimes these properties are described with the following two pairs of orthogonal complement subspaces: -* ${\mathcal C}(X)$ is the orthogonal complement of $ {\mathcal N}(X^\top )$ -* ${\mathcal R}(X)$ is the orthogonal complement ${\mathcal N}(X)$ +* ${\mathcal C}(X)$ is the orthogonal complement of ${\mathcal N}(X^\top )$ +* ${\mathcal R}(X)$ is the orthogonal complement of ${\mathcal N}(X)$ Let's do an example. @@ -281,94 +281,127 @@ import numpy as np import numpy.linalg as LA import matplotlib.pyplot as plt -rng = np.random.default_rng() +rng = np.random.default_rng(1234) ``` Having imported these modules, let's do the example. +We use a $4 \times 5$ matrix of rank $2$, so that the four subspaces have dimensions $2$, $2$, $2$, and $3$. + +Note that `np.linalg.svd` returns $V^\top$, not $V$, as its third output. + +So the row space and the null space are spanned by rows of that output. + ```{code-cell} ipython3 -np.set_printoptions(precision=2) +np.set_printoptions(precision=2, suppress=True) # Define the matrix A = np.array([[1, 2, 3, 4, 5], [2, 3, 4, 5, 6], [3, 4, 5, 6, 7], - [4, 5, 6, 7, 8], - [5, 6, 7, 8, 9]]) + [4, 5, 6, 7, 8]]) -# Compute the SVD of the matrix -U, S, V = np.linalg.svd(A,full_matrices=True) +# Compute the SVD of the matrix; the third output is V^T +U, S, VT = np.linalg.svd(A, full_matrices=True) # Compute the rank of the matrix rank = np.linalg.matrix_rank(A) -# Print the rank of the matrix print("Rank of matrix:\n", rank) print("S: \n", S) -# Compute the four fundamental subspaces -row_space = U[:, :rank] -col_space = V[:, :rank] -null_space = V[:, rank:] -left_null_space = U[:, rank:] +# Orthonormal bases for the four fundamental subspaces +col_space = U[:, :rank] # C(A), a subspace of R^4 +left_null_space = U[:, rank:] # N(A^T), a subspace of R^4 +row_space = VT[:rank, :].T # R(A), a subspace of R^5 +null_space = VT[rank:, :].T # N(A), a subspace of R^5 - -print("U:\n", U) print("Column space:\n", col_space) print("Left null space:\n", left_null_space) -print("V.T:\n", V.T) -print("Row space:\n", row_space.T) -print("Right null space:\n", null_space.T) +print("Row space:\n", row_space) +print("Null space:\n", null_space) +``` + +Let's verify that the bases have the properties that {eq}`eq:fourspaceSVD` asserts. + +```{code-cell} ipython3 +print("A @ null_space = 0: ", np.allclose(A @ null_space, 0)) +print("A.T @ left_null_space = 0: ", np.allclose(A.T @ left_null_space, 0)) +print("col_space ⟂ left_null_space:", np.allclose(col_space.T @ left_null_space, 0)) +print("row_space ⟂ null_space: ", np.allclose(row_space.T @ null_space, 0)) ``` ## Eckart-Young Theorem -Suppose that we want to construct the best rank $r$ approximation of an $m \times n$ matrix $X$. +Suppose that we want to construct the best rank $r$ approximation of an $m \times n$ matrix $X$. -By best, we mean a matrix $X_r$ of rank $r < p$ that, among all rank $r$ matrices, minimizes +By best, we mean a matrix $X_r$ of rank $r < p$ that, among all rank $r$ matrices of dimension $m \times n$, minimizes -$$ -|| X - X_r || +$$ +|| X - X_r || $$ -where $ || \cdot || $ denotes a norm of a matrix $X$ and where $X_r$ belongs to the space of all rank $r$ matrices -of dimension $m \times n$. +where $ || \cdot || $ denotes a norm of a matrix. -Three popular **matrix norms** of an $m \times n$ matrix $X$ can be expressed in terms of the singular values of $X$ +Three popular **matrix norms** of an $m \times n$ matrix $X$ can be expressed in terms of the singular values of $X$ * the **spectral** or $l^2$ norm $|| X ||_2 = \max_{||y|| \neq 0} \frac{||X y ||}{||y||} = \sigma_1$ * the **Frobenius** norm $||X ||_F = \sqrt{\sigma_1^2 + \cdots + \sigma_p^2}$ * the **nuclear** norm $ || X ||_N = \sigma_1 + \cdots + \sigma_p $ -The Eckart-Young theorem states that for each of these three norms, same rank $r$ matrix is best and that it equals +The Eckart-Young theorem states that for each of these three norms, the same rank $r$ matrix is best and that it equals $$ -\hat X_r = \sigma_1 U_1 V_1^\top + \sigma_2 U_2 V_2^\top + \cdots + \sigma_r U_r V_r^\top +\hat X_r = \sigma_1 u_1 v_1^\top + \sigma_2 u_2 v_2^\top + \cdots + \sigma_r u_r v_r^\top = U_r \Sigma_r V_r^\top $$ (eq:Ekart) -This is a very powerful theorem that says that we can take our $ m \times n $ matrix $X$ that in not full rank, and we can best approximate it by a full rank $p \times p$ matrix through the SVD. +where $U_r$ and $V_r$ consist of the first $r$ columns of $U$ and $V$, respectively, and $\Sigma_r$ is the $r \times r$ diagonal matrix of the $r$ largest singular values. + +The resulting approximation errors are + +$$ +\begin{aligned} +|| X - \hat X_r ||_2 & = \sigma_{r+1} \cr +|| X - \hat X_r ||_F & = \sqrt{\sigma_{r+1}^2 + \cdots + \sigma_p^2} \cr +|| X - \hat X_r ||_N & = \sigma_{r+1} + \cdots + \sigma_p +\end{aligned} +$$ -Moreover, if some of these $p$ singular values carry more information than others, and if we want to have the most amount of information with the least amount of data, we can take $r$ leading singular values ordered by magnitude. +This is a very powerful theorem. -We'll say more about this later when we present Principal Component Analysis. +It says that we can approximate an $m \times n$ matrix $X$ of rank $p$ by an $m \times n$ matrix of rank $r < p$ built from the $r$ largest singular values of $X$ and their singular vectors, and that the singular values that we discard measure the approximation error. + +If a few singular values are much larger than the others, a low-rank approximation captures most of $X$ while requiring us to store only $r(m + n + 1)$ numbers instead of $mn$. You can read about the Eckart-Young theorem and some of its uses [here](https://en.wikipedia.org/wiki/Low-rank_approximation). -We'll make use of this theorem when we discuss principal components analysis (PCA) and also dynamic mode decomposition (DMD). +We'll make use of this theorem when we discuss principal components analysis (PCA) below. + +It also underlies the sequel {doc}`VARs and DMDs `, where the pseudo-inverse + +$$ +\hat X_r^+ = V_r \Sigma_r^{-1} U_r^\top +$$ + +of the approximation $\hat X_r$ is a key ingredient of a dynamic mode decomposition. + +{ref}`svd_ex2` and {ref}`svd_ex3` explore these ideas. ## Full and Reduced SVD's Up to now we have described properties of a **full** SVD in which shapes of $U$, $\Sigma$, and $V$ are $\left(m, m\right)$, $\left(m, n\right)$, $\left(n, n\right)$, respectively. -There is an alternative bookkeeping convention called an **economy** or **reduced** SVD in which the shapes of $U, \Sigma$ and $V$ are different from what they are in a full SVD. +There is an alternative bookkeeping convention called an **economy** or **reduced** SVD in which the shapes of $U, \Sigma$ and $V$ are different from what they are in a full SVD. -Thus, note that because we assume that $X$ has rank $p$, there are only $p$ nonzero singular values, where $p=\textrm{rank}(X)\leq\min\left(m, n\right)$. +Because we assume that $X$ has rank $p$, there are only $p$ nonzero singular values, where $p=\textrm{rank}(X)\leq\min\left(m, n\right)$. A **reduced** SVD uses this fact to express $U$, $\Sigma$, and $V$ as matrices with shapes $\left(m, p\right)$, $\left(p, p\right)$, $\left( n, p\right)$. You can read about reduced and full SVD here +(With the option `full_matrices=False`, `numpy` returns $\min(m, n)$ columns of $U$ and $V$; when $p < \min(m,n)$, we drop the columns that correspond to zero singular values to obtain a reduced SVD in the sense used here.) + For a full SVD, $$ @@ -378,11 +411,13 @@ VV^\top & = I & \quad V^\top V = I \end{aligned} $$ -But not all these properties hold for a **reduced** SVD. +But not all these properties hold for a **reduced** SVD. -Which properties hold depend on whether we are in a **tall-skinny** case or a **short-fat** case. +For a reduced SVD, $U^\top U = I_{p \times p}$ and $V^\top V = I_{p \times p}$ always hold, but $U U^\top = I_{m \times m}$ holds only if $p = m$ and $V V^\top = I_{n \times n}$ holds only if $p = n$. - * In a **tall-skinny** case in which $m > > n$, for a **reduced** SVD +Which properties hold thus depends on whether we are in a **tall-skinny** case or a **short-fat** case. + + * In a **tall-skinny** case in which $m \gg n$ and $X$ has full column rank $p = n$, for a **reduced** SVD $$ \begin{aligned} @@ -391,19 +426,19 @@ VV^\top & = I & \quad V^\top V = I \end{aligned} $$ -* In a **short-fat** case in which $m < < n$, for a **reduced** SVD +* In a **short-fat** case in which $m \ll n$ and $X$ has full row rank $p = m$, for a **reduced** SVD $$ \begin{aligned} UU^\top & = I & \quad U^\top U = I \cr -VV^\top & = I & \quad V^\top V \neq I +VV^\top & \neq I & \quad V^\top V = I \end{aligned} $$ -When we study Dynamic Mode Decomposition below, we shall want to remember these properties when we use a reduced SVD to compute some DMD representations. +When we study dynamic mode decompositions in the sequel {doc}`VARs and DMDs `, we shall want to remember these properties because we will use a reduced SVD in a tall-skinny case. -Let's do an exercise to compare **full** and **reduced** SVD's. +Let's do an example to compare **full** and **reduced** SVD's. To review, @@ -422,13 +457,12 @@ To review, First, let's study a case in which $m = 5 > n = 2$. -(This is a small example of the **tall-skinny** case that will concern us when we study **Dynamic Mode Decompositions** below.) +(This is a small example of the **tall-skinny** case that will concern us when we study dynamic mode decompositions in the sequel.) ```{code-cell} ipython3 -import numpy as np X = rng.random((5, 2)) -U, S, V = np.linalg.svd(X,full_matrices=True) # full SVD -Uhat, Shat, Vhat = np.linalg.svd(X,full_matrices=False) # economy SVD +U, S, V = np.linalg.svd(X, full_matrices=True) # full SVD +Uhat, Shat, Vhat = np.linalg.svd(X, full_matrices=False) # economy SVD print('U, S, V =') U, S, V ``` @@ -446,8 +480,8 @@ print(f'rank of X = {rr}') **Properties:** -* Where $U$ is constructed via a full SVD, $U^\top U = I_{m\times m}$ and $U U^\top = I_{m \times m}$ -* Where $\hat U$ is constructed via a reduced SVD, although $\hat U^\top \hat U = I_{p\times p}$, it happens that $\hat U \hat U^\top \neq I_{m \times m}$ +* Where $U$ is constructed via a full SVD, $U^\top U = I_{m\times m}$ and $U U^\top = I_{m \times m}$ +* Where $\hat U$ is constructed via a reduced SVD, although $\hat U^\top \hat U = I_{p\times p}$, it happens that $\hat U \hat U^\top \neq I_{m \times m}$ We illustrate these properties for our example with the following code cells. @@ -467,30 +501,25 @@ UhatUhatT, UhatTUhat ``` - - **Remarks:** -The cells above illustrate the application of the `full_matrices=True` and `full_matrices=False` options. -Using `full_matrices=False` returns a reduced singular value decomposition. - -The **full** and **reduced** SVD's both accurately decompose an $m \times n$ matrix $X$ - -When we study Dynamic Mode Decompositions below, it will be important for us to remember the preceding properties of full and reduced SVD's in such tall-skinny cases. +The cells above illustrate the application of the `full_matrices=True` and `full_matrices=False` options. +Using `full_matrices=False` returns a reduced singular value decomposition. +The **full** and **reduced** SVD's both accurately decompose an $m \times n$ matrix $X$. +When we study dynamic mode decompositions in the sequel, it will be important for us to remember the preceding properties of full and reduced SVD's in such tall-skinny cases. Now let's turn to a short-fat case. -To illustrate this case, we'll set $m = 2 < 5 = n $ and compute both full and reduced SVD's. +To illustrate this case, we'll set $m = 2 < 5 = n$ and compute both full and reduced SVD's. ```{code-cell} ipython3 -import numpy as np X = rng.random((2, 5)) -U, S, V = np.linalg.svd(X,full_matrices=True) # full SVD -Uhat, Shat, Vhat = np.linalg.svd(X,full_matrices=False) # economy SVD +U, S, V = np.linalg.svd(X, full_matrices=True) # full SVD +Uhat, Shat, Vhat = np.linalg.svd(X, full_matrices=False) # economy SVD print('U, S, V = ') U, S, V ``` @@ -499,6 +528,7 @@ U, S, V print('Uhat, Shat, Vhat = ') Uhat, Shat, Vhat ``` + Let's verify that our reduced SVD accurately represents $X$ ```{code-cell} ipython3 @@ -523,66 +553,72 @@ Q & = U V^\top \end{aligned} $$ -Here - -* $S$ is an $m \times m$ **symmetric** matrix -* $Q$ is an $m \times n$ **orthogonal** matrix - and in our reduced SVD * $U$ is an $m \times p$ orthonormal matrix * $\Sigma$ is a $p \times p$ diagonal matrix -* $V$ is an $n \times p$ orthonormal +* $V$ is an $n \times p$ orthonormal matrix + +Because $U^\top U = I$, we have $SQ = U \Sigma U^\top U V^\top = U \Sigma V^\top = X$. + +Here + +* $S$ is an $m \times m$ **symmetric positive semidefinite** matrix +* $Q$ is an $m \times n$ matrix with $Q Q^\top = U U^\top$; so when $p = m$, as for a short-fat $X$ with full row rank, $Q Q^\top = I$ and $Q$ has orthonormal rows ## Application: Principal Components Analysis (PCA) -Let's begin with a case in which $n >> m$, so that we have many more individuals $n$ than attributes $m$. +Let's begin with a case in which $n \gg m$, so that we have many more individuals $n$ than attributes $m$. -The matrix $X$ is **short and fat** in an $n >> m$ case as opposed to a **tall and skinny** case with $m > > n $ to be discussed later. +The matrix $X$ is **short and fat** in an $n \gg m$ case as opposed to a **tall and skinny** case with $m \gg n$ to be discussed in the sequel. -We regard $X$ as an $m \times n$ matrix of **data**: +We regard $X$ as an $m \times n$ matrix of **data**: $$ X = \begin{bmatrix} X_1 \mid X_2 \mid \cdots \mid X_n\end{bmatrix} $$ -where for $j = 1, \ldots, n$ the column vector $X_j = \begin{bmatrix}x_{1j}\\x_{2j}\\\vdots\\x_{mj}\end{bmatrix}$ is a vector of observations on variables $\begin{bmatrix}X_1\\X_2\\\vdots\\X_m\end{bmatrix}$. +where for $j = 1, \ldots, n$ the column vector $X_j = \begin{bmatrix}x_{1j}\\x_{2j}\\\vdots\\x_{mj}\end{bmatrix}$ is a vector of observations on variables $1, 2, \ldots, m$. In a **time series** setting, we would think of columns $j$ as indexing different __times__ at which random variables are observed, while rows index different random variables. -In a **cross-section** setting, we would think of columns $j$ as indexing different __individuals__ for which random variables are observed, while rows index different **attributes**. +In a **cross-section** setting, we would think of columns $j$ as indexing different __individuals__ for which random variables are observed, while rows index different **attributes**. + +As we have seen before, the SVD is a way to decompose a matrix into useful components, just like polar decomposition, eigendecomposition, and many others. -As we have seen before, the SVD is a way to decompose a matrix into useful components, just like polar decomposition, eigendecomposition, and many others. +PCA, on the other hand, is a method that builds on the SVD to analyze data. -PCA, on the other hand, is a method that builds on the SVD to analyze data. The goal is to apply certain steps, to help better visualize patterns in data, using statistical tools to capture the most important patterns in data. +Its goal is to find a few linear combinations of the variables that capture the most important patterns of variation in the data. -**Step 1: Standardize the data:** +**Step 1: Center the data:** -Because our data matrix may hold variables of different units and scales, we first need to standardize the data. +Because we are interested in variation of the data around their means, we first subtract sample means. -First by computing the average of each row of $X$. +(Because our data matrix may hold variables of different units and scales, we might also **standardize** the data by dividing each row by its sample standard deviation.) + +We first compute the average of each row of $X$ $$ \bar{X_i}= \frac{1}{n} \sum_{j = 1}^{n} x_{ij} $$ -We then create an average matrix out of these means: +We then create a matrix of these means: $$ -\bar{X} = \begin{bmatrix} \bar{X_1} \\ \bar{X_2} \\ \ldots \\ \bar{X_m}\end{bmatrix}\begin{bmatrix}1 \mid 1 \mid \cdots \mid 1 \end{bmatrix} +\bar{X} = \begin{bmatrix} \bar{X_1} \\ \bar{X_2} \\ \vdots \\ \bar{X_m}\end{bmatrix}\begin{bmatrix}1 \mid 1 \mid \cdots \mid 1 \end{bmatrix} $$ -And subtract out of the original matrix to create a mean centered matrix: +and subtract it from the original matrix to create a mean-centered matrix: $$ B = X - \bar{X} $$ -**Step 2: Compute the covariance matrix:** +**Step 2: Compute the covariance matrix:** -Then because we want to extract the relationships between variables rather than just their magnitude, in other words, we want to know how they can explain each other, we compute the covariance matrix of $B$. +Because we want to extract relationships between variables rather than just their magnitudes -- in other words, we want to know how they can explain each other -- we compute the sample covariance matrix of $B$. $$ C = \frac{1}{n} BB^{\top} @@ -590,7 +626,7 @@ $$ **Step 3: Decompose the covariance matrix and arrange the singular values:** -Since the matrix $C$ is positive definite, we can eigendecompose it, find its eigenvalues, and rearrange the eigenvalue and eigenvector matrices in a decreasing order. +Since the matrix $C$ is symmetric and positive semidefinite, we can eigendecompose it, find its eigenvalues, and arrange the eigenvalue and eigenvector matrices in decreasing order. The eigendecomposition of $C$ can be found by decomposing $B$ instead. Since $B$ is not a square matrix, we obtain an SVD of $B$: @@ -602,39 +638,51 @@ B B^\top &= U \Sigma V^\top (U \Sigma V^{\top})^{\top}\\ \end{aligned} $$ +so that + $$ C = \frac{1}{n} U \Sigma \Sigma^\top U^\top $$ -We can then rearrange the columns in the matrices $U$ and $\Sigma$ so that the singular values are in decreasing order. +Singular values are conventionally arranged in decreasing order, as `numpy` does. **Step 4: Select singular values, (optional) truncate the rest:** -We can now decide how many singular values to pick, based on how much variance you want to retain. (e.g., retaining 95% of the total variance). +We can now decide how many singular values to keep, based on how much variance we want to retain (e.g., retaining 95% of the total variance). -We can obtain the percentage by calculating the variance contained in the leading $r$ factors divided by the variance in total: +We can obtain the percentage by calculating the variance contained in the leading $r$ factors divided by the total variance: $$ \frac{\sum_{i = 1}^{r} \sigma^2_{i}}{\sum_{i = 1}^{p} \sigma^2_{i}} $$ +By the Eckart-Young theorem, one minus this ratio is also the squared relative Frobenius-norm error $||B - \hat B_r||_F^2 / ||B||_F^2$ of the best rank $r$ approximation $\hat B_r$ of $B$. + **Step 5: Create the Score Matrix:** +The matrix of **principal components** (or **scores**) is + $$ \begin{aligned} -T&= BV \cr -&= U\Sigma V^\top V \cr -&= U\Sigma +T & = U^\top B \cr +& = U^\top U \Sigma V^\top \cr +& = \Sigma V^\top \end{aligned} $$ +Its $k$th row is the sequence, across observations $j = 1, \ldots, n$, of the $k$th principal component. + +Keeping only the first $r$ rows gives the $r \times n$ matrix $T_r = U_r^\top B$ of the first $r$ principal components, and $U_r T_r = U_r \Sigma_r V_r^\top = \hat B_r$ is the Eckart-Young approximation of $B$. + +In the sequel {doc}`VARs and DMDs `, the vector $U_r^\top X_t$ of the first $r$ principal components of a time $t$ cross section serves as the regressor in a **principal components regression** that forecasts next period's cross section. + ## Relationship of PCA to SVD -To relate an SVD to a PCA of data set $X$, first construct the SVD of the data matrix $X$: +To relate an SVD to a PCA of data set $X$, first construct the SVD of the data matrix $X$. -Let’s assume that sample means of all variables are zero, so we don't need to standardize our matrix. +Let’s assume that sample means of all variables are zero, so we don't need to center our matrix. $$ X = U \Sigma V^\top = \sigma_1 U_1 V_1^\top + \sigma_2 U_2 V_2^\top + \cdots + \sigma_p U_p V_p^\top @@ -647,7 +695,7 @@ U=\begin{bmatrix}U_1|U_2|\ldots|U_m\end{bmatrix} $$ $$ -V^\top = \begin{bmatrix}V_1^\top \\V_2^\top \\\ldots\\V_n^\top \end{bmatrix} +V^\top = \begin{bmatrix}V_1^\top \\V_2^\top \\\vdots\\V_n^\top \end{bmatrix} $$ In equation {eq}`eq:PCA1`, each of the $m \times n$ matrices $U_{j}V_{j}^\top $ is evidently @@ -659,27 +707,27 @@ $$ X = \sigma_1 \begin{bmatrix}U_{11}V_{1}^\top \\U_{21}V_{1}^\top \\\cdots\\U_{m1}V_{1}^\top \\\end{bmatrix} + \sigma_2\begin{bmatrix}U_{12}V_{2}^\top \\U_{22}V_{2}^\top \\\cdots\\U_{m2}V_{2}^\top \\\end{bmatrix}+\ldots + \sigma_p\begin{bmatrix}U_{1p}V_{p}^\top \\U_{2p}V_{p}^\top \\\cdots\\U_{mp}V_{p}^\top \\\end{bmatrix} $$ (eq:PCA2) -Here is how we would interpret the objects in the matrix equation {eq}`eq:PCA2` in +Here is how we would interpret the objects in the matrix equation {eq}`eq:PCA2` in a time series context: -* $ \textrm{for each} \ k=1, \ldots, n $, the object $\lbrace V_{kj} \rbrace_{j=1}^n$ is a time series for the $k$th **principal component** +* for each $k=1, \ldots, p$, the object $\lbrace V_{jk} \rbrace_{j=1}^n$ (the $k$th column of $V$) is a time series for the $k$th **principal component**, normalized to have unit length; the unnormalized principal component is $\sigma_k V_k^\top = U_k^\top X$ -* $U_j = \begin{bmatrix}U_{1k}\\U_{2k}\\\ldots\\U_{mk}\end{bmatrix} \ k=1, \ldots, m$ -is a vector of **loadings** of variables $X_i$ on the $k$th principal component, $i=1, \ldots, m$ +* $U_k = \begin{bmatrix}U_{1k}\\U_{2k}\\\vdots\\U_{mk}\end{bmatrix}, \ k=1, \ldots, p$, +is a vector of **loadings** of variables $X_i$ on the $k$th principal component, $i=1, \ldots, m$ -* $\sigma_k $ for each $k=1, \ldots, p$ is the strength of $k$th **principal component**, where strength means contribution to the overall covariance of $X$. +* $\sigma_k $ for each $k=1, \ldots, p$ is the strength of the $k$th **principal component**, where strength means contribution to the overall covariance of $X$: the $k$th principal component contributes $\sigma_k^2 / n$ to the sum of the sample variances of the $m$ variables ## PCA with Eigenvalues and Eigenvectors -We now use an eigen decomposition of a sample covariance matrix to do PCA. +We now use an eigen decomposition of a sample covariance matrix to do PCA. Let $X_{m \times n}$ be our $m \times n$ data matrix. Let's assume that sample means of all variables are zero. -We can assure this by **pre-processing** the data by subtracting sample means. +We can assure this by **pre-processing** the data by subtracting sample means. -Define a sample covariance matrix $\Omega$ as +Define a matrix $\Omega$ that is proportional to the sample covariance matrix (it equals $n$ times it) as $$ \Omega = XX^\top @@ -693,7 +741,7 @@ $$ Here -* $P$ is $m×m$ matrix of eigenvectors of $\Omega$ +* $P$ is an $m×m$ orthogonal matrix of eigenvectors of $\Omega$ * $\Lambda$ is a diagonal matrix of eigenvalues of $\Omega$ @@ -706,7 +754,7 @@ $$ where $$ -\epsilon = P^{-1} X +\epsilon = P^{-1} X = P^\top X $$ and @@ -721,44 +769,47 @@ $$ XX^\top =P\Lambda P^\top . $$ (eq:XXo) -It follows that we can represent the data matrix $X$ as +It follows that we can represent the data matrix $X$ as -\begin{equation*} -X=\begin{bmatrix}X_1|X_2|\ldots|X_m\end{bmatrix} =\begin{bmatrix}P_1|P_2|\ldots|P_m\end{bmatrix} -\begin{bmatrix}\epsilon_1\\\epsilon_2\\\ldots\\\epsilon_m\end{bmatrix} +$$ +X =\begin{bmatrix}P_1|P_2|\ldots|P_m\end{bmatrix} +\begin{bmatrix}\epsilon_1\\\epsilon_2\\\vdots\\\epsilon_m\end{bmatrix} = P_1\epsilon_1+P_2\epsilon_2+\ldots+P_m\epsilon_m -\end{equation*} +$$ + +where $P_j$ is the $j$th column of $P$ and $\epsilon_j$ is the $1 \times n$ $j$th row of $\epsilon$. +To reconcile the preceding representation with the PCA that we had obtained earlier through the SVD, we first note that $\epsilon_j \epsilon_j^\top = \lambda_j \equiv \sigma^2_j$. -To reconcile the preceding representation with the PCA that we had obtained earlier through the SVD, we first note that $\epsilon_j^2=\lambda_j\equiv\sigma^2_j$. +Now for each $j$ with $\lambda_j > 0$ define $\tilde{\epsilon}_j = \frac{\epsilon_j}{\sqrt{\lambda_j}}$, +which implies that $\tilde{\epsilon}_j\tilde{\epsilon}_j^\top =1$. -Now define $\tilde{\epsilon_j} = \frac{\epsilon_j}{\sqrt{\lambda_j}}$, -which implies that $\tilde{\epsilon}_j\tilde{\epsilon}_j^\top =1$. +Terms with $\lambda_j = 0$ vanish because then $\epsilon_j = 0$. Therefore $$ \begin{aligned} -X&=\sqrt{\lambda_1}P_1\tilde{\epsilon_1}+\sqrt{\lambda_2}P_2\tilde{\epsilon_2}+\ldots+\sqrt{\lambda_m}P_m\tilde{\epsilon_m}\\ -&=\sigma_1P_1\tilde{\epsilon_2}+\sigma_2P_2\tilde{\epsilon_2}+\ldots+\sigma_mP_m\tilde{\epsilon_m} , +X&=\sqrt{\lambda_1}P_1\tilde{\epsilon}_1+\sqrt{\lambda_2}P_2\tilde{\epsilon}_2+\ldots+\sqrt{\lambda_p}P_p\tilde{\epsilon}_p\\ +&=\sigma_1P_1\tilde{\epsilon}_1+\sigma_2P_2\tilde{\epsilon}_2+\ldots+\sigma_pP_p\tilde{\epsilon}_p , \end{aligned} $$ -which agrees with +which agrees with $$ -X=\sigma_1U_1{V_1}^{T}+\sigma_2 U_2{V_2}^{T}+\ldots+\sigma_{r} U_{r}{V_{r}}^{T} +X=\sigma_1U_1{V_1}^{\top}+\sigma_2 U_2{V_2}^{\top}+\ldots+\sigma_{p} U_{p}{V_{p}}^{\top} $$ -provided that we set +provided that we set -* $U_j=P_j$ (a vector of loadings of variables on principal component $j$) +* $U_j=P_j$ (a vector of loadings of variables on principal component $j$) -* ${V_k}^{T}=\tilde{\epsilon_k}$ (the $k$th principal component) +* ${V_k}^{\top}=\tilde{\epsilon}_k$ (the $k$th principal component) -Because there are alternative algorithms for computing $P$ and $U$ for given a data matrix $X$, depending on algorithms used, we might have sign differences or different orders of eigenvectors. +Because there are alternative algorithms for computing $P$ and $U$ for a given data matrix $X$, depending on the algorithms used, we might have sign differences or different orders of eigenvectors. -We can resolve such ambiguities about $U$ and $P$ by +We can resolve such ambiguities about $U$ and $P$ by 1. sorting eigenvalues and singular values in descending order 2. imposing positive diagonals on $P$ and $U$ and adjusting signs in $V^\top $ accordingly @@ -767,7 +818,7 @@ We can resolve such ambiguities about $U$ and $P$ by To pull things together, it is useful to assemble and compare some formulas presented above. -First, consider an SVD of an $m \times n$ matrix: +First, consider an SVD of an $m \times n$ matrix: $$ X = U\Sigma V^\top @@ -778,14 +829,14 @@ Compute: $$ \begin{aligned} XX^\top &=U\Sigma V^\top V\Sigma^\top U^\top \cr -&\equiv U\Sigma\Sigma^\top U^\top \cr +&= U\Sigma\Sigma^\top U^\top \cr &\equiv U\Lambda U^\top \end{aligned} $$ (eq:XXcompare) Compare representation {eq}`eq:XXcompare` with equation {eq}`eq:XXo` above. -Evidently, $U$ in the SVD is the matrix $P$ of +Evidently, $U$ in the SVD is the matrix $P$ of eigenvectors of $XX^\top $ and $\Sigma \Sigma^\top $ is the matrix $\Lambda$ of eigenvalues. Second, let's compute @@ -799,7 +850,7 @@ $$ -Thus, the matrix $V$ in the SVD is the matrix of eigenvectors of $X^\top X$ +Thus, the matrix $V$ in the SVD is the matrix of eigenvectors of $X^\top X$. Summarizing and fitting things together, we have the eigen decomposition of the sample covariance matrix @@ -818,7 +869,7 @@ $$ where $U$ is an orthogonal matrix. -Thus, $P = U$ and we have the representation of $X$ +Thus, $P = U$ (up to the ordering and signs of columns) and we have the representation of $X$ $$ X = P \epsilon = U \Sigma V^\top @@ -830,6 +881,8 @@ $$ U^\top X = \Sigma V^\top = \epsilon $$ +so that $\epsilon$ is the matrix of principal components (scores) that we constructed in Step 5 above. + Note that the preceding implies that $$ @@ -838,7 +891,7 @@ $$ so that everything fits together. -Below we define a class `DecomAnalysis` that wraps PCA and SVD for a given a data matrix `X`. +Below we define a class `DecomAnalysis` that wraps PCA and SVD for a given data matrix `X`. ```{code-cell} ipython3 class DecomAnalysis: @@ -860,7 +913,7 @@ class DecomAnalysis: if r_component: self.r_component = r_component else: - self.r_component = self.m + self.r_component = min(self.m, self.n) def pca(self): @@ -877,7 +930,7 @@ class DecomAnalysis: self.explained_ratio_pca = np.cumsum(self.𝜆) / self.𝜆.sum() - # compute the N by T matrix of principal components + # compute the m by n matrix of principal components self.𝜖 = self.P.T @ self.X P = self.P[:, :self.r_component] @@ -941,7 +994,7 @@ def diag_sign(A): return D ``` -We also define a function that prints out information so that we can compare decompositions +We also define a function that prints out information so that we can compare decompositions obtained by different algorithms. ```{code-cell} ipython3 @@ -958,61 +1011,97 @@ def compare_pca_svd(da): print(f'σ^2 = {da.σ**2}\n') print('\n') + k = da.r + # loading matrices fig, axs = plt.subplots(1, 2, figsize=(14, 5)) plt.suptitle('loadings') - axs[0].plot(da.P.T) + axs[0].plot(da.P[:, :k]) axs[0].set_title('P') - axs[0].set_xlabel('m') - axs[1].plot(da.U.T) + axs[0].set_xlabel('variable $i$') + axs[1].plot(da.U[:, :k]) axs[1].set_title('U') - axs[1].set_xlabel('m') + axs[1].set_xlabel('variable $i$') plt.show() # principal components fig, axs = plt.subplots(1, 2, figsize=(14, 5)) plt.suptitle('principal components') - axs[0].plot(da.ε.T) + axs[0].plot(da.ε[:k, :].T) axs[0].set_title('ε') axs[0].set_xlabel('n') - axs[1].plot(da.VT[:da.r, :].T * np.sqrt(da.λ)) + axs[1].plot(da.VT[:k, :].T * np.sqrt(da.λ[:k])) axs[1].set_title(r'$V^\top *\sqrt{\lambda}$') axs[1].set_xlabel('n') plt.show() ``` +Let's apply these tools to a short-fat data matrix with $m = 5$ variables and $n = 200$ observations that is generated by two factors plus noise. + +```{code-cell} ipython3 +m, n = 5, 200 +G = rng.normal(size=(m, 2)) # loadings on two factors +Z = rng.normal(size=(2, n)) # factors +X = G @ Z + 0.2 * rng.normal(size=(m, n)) +X = X - X.mean(axis=1, keepdims=True) # subtract sample means + +da = DecomAnalysis(X) +compare_pca_svd(da) +``` + +The two algorithms produce the same loadings and the same principal components. + +The cumulative fractions of variance explained reveal that two principal components account for almost all of the variation in $X$, which is consistent with how we generated the data. + +```{code-cell} ipython3 +print("explained ratio (PCA):", da.explained_ratio_pca) +print("explained ratio (SVD):", da.explained_ratio_svd) +``` + +For an example of PCA applied to analyzing the structure of intelligence tests see the lecture {doc}`Multivariate Normal Distribution `. + +Look at parts of that lecture that describe and illustrate the classic factor analysis model. + +As mentioned earlier, in a sequel to this lecture, {doc}`VARs and DMDs `, we'll describe how SVD's provide ways rapidly to compute reduced-order approximations to first-order Vector Autoregressions (VARs). + ## Exercises ```{exercise} :label: svd_ex1 -In Ordinary Least Squares (OLS), we learn to compute $ \hat{\beta} = (X^\top X)^{-1} X^\top y $, but there are cases such as when we have colinearity or an underdetermined system: **short fat** matrix. +In ordinary least squares (OLS), we learn to compute $ \hat{\beta} = (X^\top X)^{-1} X^\top y $, where the matrix $X$ of regressors has one row for each observation and one column for each regressor. -In these cases, the $ (X^\top X) $ matrix is not not invertible (its determinant is zero) or ill-conditioned (its determinant is very close to zero). +But the matrix $X^\top X$ is not invertible (its determinant is zero) when regressors are perfectly collinear or when there are more regressors than observations, so that $X$ is **short and fat** and the system $X \beta = y$ is underdetermined. -What we can do instead is to create what is called a [pseudoinverse](https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse), a full rank approximation of the inverted matrix so we can compute $ \hat{\beta} $ with it. +And $X^\top X$ is **ill-conditioned** when regressors are nearly collinear, meaning that the ratio of its largest to its smallest eigenvalue is very large. -Thinking in terms of the Eckart-Young theorem, build the pseudoinverse matrix $ X^{+} $ and use it to compute $ \hat{\beta} $. +What we can do instead is to use what is called a [pseudoinverse](https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse) $X^+$ of $X$ in place of $(X^\top X)^{-1} X^\top$. +1. Use an SVD of $X$ to build the pseudoinverse matrix $X^{+}$ and use it to compute $\hat{\beta} = X^+ y$. +1. Show that when $X$ has linearly independent columns (an **over-determined** system), $X^+ = (X^\top X)^{-1} X^\top$, so that $\hat \beta$ is the usual OLS estimator. +1. Show that when $X$ has linearly independent rows (an **under-determined** system), $\hat \beta$ solves $X \beta = y$ exactly and has the smallest norm among all solutions. +1. Verify these claims numerically for a $5 \times 8$ matrix $X$ of random numbers. ``` ```{solution-start} svd_ex1 :class: dropdown ``` -We can use SVD to compute the pseudoinverse: +We can use an SVD to compute the pseudoinverse. + +Starting from a reduced SVD $$ X = U \Sigma V^\top $$ -inverting $X$, we have: +in which $\Sigma$ is the $p \times p$ diagonal matrix of positive singular values, the pseudoinverse is $$ -X^{+} = V \Sigma^{+} U^\top +X^{+} = V \Sigma^{-1} U^\top $$ -where: +Equivalently, in terms of a full SVD, $X^+ = V \Sigma^+ U^\top$, where the $n \times m$ matrix $\Sigma^+$ is formed by transposing $\Sigma$ and replacing each positive singular value $\sigma_j$ by $1/\sigma_j$: $$ \Sigma^{+} = \begin{bmatrix} @@ -1027,15 +1116,201 @@ $$ and finally: $$ -\hat{\beta} = X^{+}y = V \Sigma^{+} U^\top y +\hat{\beta} = X^{+}y = V \Sigma^{-1} U^\top y +$$ + +When $X$ has linearly independent columns, $V$ is a square orthogonal matrix, so + +$$ +(X^\top X)^{-1} X^\top = (V \Sigma^2 V^\top)^{-1} V \Sigma U^\top = V \Sigma^{-2} V^\top V \Sigma U^\top = V \Sigma^{-1} U^\top = X^+ . +$$ + +When $X$ has linearly independent rows, $U$ is a square orthogonal matrix, so + +$$ +X \hat \beta = U \Sigma V^\top V \Sigma^{-1} U^\top y = U U^\top y = y . $$ +Every other solution of $X \beta = y$ has the form $\beta = \hat \beta + \eta$ with $X \eta = 0$. + +The vector $\hat \beta = V (\Sigma^{-1} U^\top y)$ lies in the row space of $X$, which is spanned by the columns of $V$, while $\eta$ lies in the null space of $X$. + +These two subspaces are orthogonal, so $\|\beta\|^2 = \|\hat \beta\|^2 + \|\eta\|^2 \geq \|\hat \beta\|^2$. + +Here is a numerical check. + +```{code-cell} ipython3 +n_obs, k = 5, 8 +X = rng.normal(size=(n_obs, k)) +y = rng.normal(size=n_obs) + +U, σ, VT = np.linalg.svd(X, full_matrices=False) +X_plus = VT.T @ np.diag(1 / σ) @ U.T +β_hat = X_plus @ y + +print("X_plus equals np.linalg.pinv(X):", np.allclose(X_plus, np.linalg.pinv(X))) +print("X β_hat = y: ", np.allclose(X @ β_hat, y)) + +# another exact solution: add an element of the null space of X +β_other = β_hat + (np.eye(k) - X_plus @ X) @ rng.normal(size=k) +print("X β_other = y: ", np.allclose(X @ β_other, y)) +print(f"norms: |β_hat| = {LA.norm(β_hat):.3f}, |β_other| = {LA.norm(β_other):.3f}") + +# over-determined case: pseudo-inverse reproduces OLS +X_tall = rng.normal(size=(8, 5)) +print("pinv equals (X^T X)^{-1} X^T: ", + np.allclose(np.linalg.pinv(X_tall), LA.inv(X_tall.T @ X_tall) @ X_tall.T)) +``` + +```{solution-end} +``` + +```{exercise} +:label: svd_ex2 + +This exercise verifies the Eckart-Young theorem numerically. + +Construct a $60 \times 40$ matrix $X$ equal to a rank-$4$ matrix plus noise, for example + +`X = rng.normal(size=(60, 4)) @ rng.normal(size=(4, 40)) + 0.3 * rng.normal(size=(60, 40))` + +For $r = 1, 2, 4, 8$ + +1. compute the Eckart-Young approximation $\hat X_r$ and verify the formulas for $|| X - \hat X_r ||$ in the spectral, Frobenius, and nuclear norms +1. compare $|| X - \hat X_r ||_F$ with $|| X - Q Q^\top X ||_F$ for $500$ matrices $Q$ with $r$ orthonormal columns that span randomly drawn subspaces (for each $Q$, $Q Q^\top X$ is the best approximation of $X$ whose columns lie in the column space of $Q$) + +Finally, plot the fraction $\sum_{i \leq r} \sigma_i^2 / \sum_i \sigma_i^2$ against $r$ and relate it to Step 4 of the PCA recipe above. +``` + +```{solution-start} svd_ex2 +:class: dropdown +``` + +Here is one solution. + +```{code-cell} ipython3 +m, n = 60, 40 +X = rng.normal(size=(m, 4)) @ rng.normal(size=(4, n)) + 0.3 * rng.normal(size=(m, n)) +U, σ, VT = np.linalg.svd(X, full_matrices=False) + +print(" r spectral Frobenius nuclear |X - X_r|_F best random |X - QQ'X|_F") +for r in (1, 2, 4, 8): + X_r = U[:, :r] @ np.diag(σ[:r]) @ VT[:r, :] + E = X - X_r + checks = (np.isclose(LA.norm(E, 2), σ[r]), + np.isclose(LA.norm(E, 'fro'), np.sqrt(np.sum(σ[r:]**2))), + np.isclose(LA.norm(E, 'nuc'), np.sum(σ[r:]))) + best_random = np.inf + for _ in range(500): + Q = LA.qr(rng.normal(size=(m, r)))[0] + best_random = min(best_random, LA.norm(X - Q @ Q.T @ X, 'fro')) + print(f"{r:2d} {str(checks[0]):8s} {str(checks[1]):9s} {str(checks[2]):8s}" + f" {LA.norm(E, 'fro'):10.2f} {best_random:10.2f}") +``` + +All three error formulas hold, and no randomly drawn rank-$r$ approximation comes close to the Eckart-Young approximation. + +```{code-cell} ipython3 +share = np.cumsum(σ**2) / np.sum(σ**2) + +fig, ax = plt.subplots() +ax.plot(np.arange(1, len(σ) + 1), share, 'o-') +ax.set_xlabel('$r$') +ax.set_ylabel('fraction of $||X||_F^2$ captured by $\\hat X_r$') +plt.show() +``` + +The fraction rises steeply until $r = 4$, the rank of the signal part of $X$, and then flattens out. + +This fraction is exactly the "fraction of variance retained" that Step 4 of the PCA recipe uses to choose how many principal components to keep, and one minus it is the squared relative Frobenius error of $\hat X_r$. + ```{solution-end} ``` +```{exercise} +:label: svd_ex3 + +This exercise builds a bridge to the sequel {doc}`VARs and DMDs `. + +Simulate a **tall and skinny** time series data set $\tilde X$ with $m = 100$ variables and $61$ periods in which + +$$ +\tilde X_t = G z_t + 0.5 v_t, \qquad z_{t+1} = \begin{bmatrix} 0.9 & 0 \cr 0 & 0.6 \end{bmatrix} z_t + w_{t+1} +$$ + +where $G$ is a $100 \times 2$ matrix of standard normal random numbers and $v_t$ and $w_{t+1}$ are standard normal random vectors. + +Subtract the sample mean of each variable, then use the first $41$ periods to form $X = \begin{bmatrix} \tilde X_1 \mid \cdots \mid \tilde X_{40} \end{bmatrix}$ and $X' = \begin{bmatrix} \tilde X_2 \mid \cdots \mid \tilde X_{41} \end{bmatrix}$. + +1. Compute the least squares estimator $\hat A = X' X^+$ of the matrix $A$ in the first-order VAR $\tilde X_{t+1} = A \tilde X_t + \epsilon_{t+1}$ and verify that it fits perfectly: $\hat A X = X'$. +1. Let $\hat X_2$ be the Eckart-Young rank-$2$ approximation of $X$. Verify numerically that $\hat X_2^+ = V_2 \Sigma_2^{-1} U_2^\top$ and compute the rank-$2$ estimator $\hat A_2 = X' \hat X_2^+$. +1. Let $F = U_2^\top X$ be the $2 \times 40$ matrix of the first two principal components of $X$. Verify that $\hat A_2 = X' F^\top (F F^\top)^{-1} U_2^\top$, i.e., that $\hat A_2 \tilde X_t$ equals the fitted value from a regression of $\tilde X_{t+1}$ on the first two principal components of $\tilde X_t$. +1. Use the remaining $20$ periods to compare one-step-ahead out-of-sample forecast errors $\|X'_{\rm out} - \check A X_{\rm out}\|_F / \|X'_{\rm out}\|_F$ for $\check A = \hat A$ and $\check A = \hat A_2$. Repeat for a few random seeds. +``` + +```{solution-start} svd_ex3 +:class: dropdown +``` + +Here is one solution. + +```{code-cell} ipython3 +def simulate_tall(m=100, T=61, σ_v=0.5, seed=1): + rng = np.random.default_rng(seed) + G = rng.normal(size=(m, 2)) + a = np.array([0.9, 0.6]) + z = np.zeros((2, T)) + for t in range(T - 1): + z[:, t+1] = a * z[:, t] + rng.normal(size=2) + X_tilde = G @ z + σ_v * rng.normal(size=(m, T)) + return X_tilde - X_tilde.mean(axis=1, keepdims=True) + +def rel_error(A, X, X_prime): + return LA.norm(X_prime - A @ X) / LA.norm(X_prime) + +X_tilde = simulate_tall() +X, X_prime = X_tilde[:, :40], X_tilde[:, 1:41] # estimation sample +X_out, X_prime_out = X_tilde[:, 40:-1], X_tilde[:, 41:] # hold-out sample + +# 1. minimum-norm least squares estimator +A_hat = X_prime @ np.linalg.pinv(X) +print("A_hat X = X':", np.allclose(A_hat @ X, X_prime)) + +# 2. pseudo-inverse of the Eckart-Young approximation +U, σ, VT = np.linalg.svd(X, full_matrices=False) +r = 2 +U_r, σ_r, V_r = U[:, :r], σ[:r], VT[:r, :].T +X_r = U_r @ np.diag(σ_r) @ V_r.T +X_r_plus = V_r @ np.diag(1 / σ_r) @ U_r.T +print("pinv(X_r) = V_r Σ_r^{-1} U_r^T:", np.allclose(np.linalg.pinv(X_r, rcond=1e-10), X_r_plus)) +A_r = X_prime @ X_r_plus + +# 3. principal components regression +F = U_r.T @ X +A_pcr = X_prime @ F.T @ LA.inv(F @ F.T) @ U_r.T +print("A_r equals principal components regression:", np.allclose(A_r, A_pcr)) + +# 4. in-sample and out-of-sample fit +print(f"\nin-sample: A_hat {rel_error(A_hat, X, X_prime):.3f}, A_r {rel_error(A_r, X, X_prime):.3f}") +print(f"out-of-sample: A_hat {rel_error(A_hat, X_out, X_prime_out):.3f}, A_r {rel_error(A_r, X_out, X_prime_out):.3f}") + +print("\nout-of-sample errors for other seeds") +for seed in range(2, 7): + X_tilde = simulate_tall(seed=seed) + X, X_prime = X_tilde[:, :40], X_tilde[:, 1:41] + X_out, X_prime_out = X_tilde[:, 40:-1], X_tilde[:, 41:] + U, σ, VT = np.linalg.svd(X, full_matrices=False) + A_r = X_prime @ VT[:r, :].T @ np.diag(1 / σ[:r]) @ U[:, :r].T + A_hat = X_prime @ np.linalg.pinv(X) + print(f"seed {seed}: A_hat {rel_error(A_hat, X_out, X_prime_out):.3f}," + f" A_r {rel_error(A_r, X_out, X_prime_out):.3f}") +``` + +The minimum-norm estimator $\hat A$ fits the estimation sample perfectly but forecasts poorly out of sample, because the pseudo-inverse $X^+ = V \Sigma^{-1} U^\top$ puts large weights $1/\sigma_j$ on directions with small singular values that mostly contain noise. -For an example PCA applied to analyzing the structure of intelligence tests see this lecture {doc}`Multivariable Normal Distribution `. +The rank-$2$ estimator $\hat A_2$ discards those directions and forecasts better out of sample. -Look at parts of that lecture that describe and illustrate the classic factor analysis model. +The estimator $\hat A_2 = X' \hat X_2^+$ is exactly the estimator that a **dynamic mode decomposition** computes; the sequel {doc}`VARs and DMDs ` studies it in detail. -As mentioned earlier, in a sequel to this lecture about {doc}`Dynamic Mode Decompositions `, we'll describe how SVD's provide ways rapidly to compute reduced-order approximations to first-order Vector Autoregressions (VARs). +```{solution-end} +``` diff --git a/lectures/var_dmd.md b/lectures/var_dmd.md index ee8e8c2bd..67196c145 100644 --- a/lectures/var_dmd.md +++ b/lectures/var_dmd.md @@ -13,54 +13,72 @@ kernelspec: # VARs and DMDs -This lecture applies computational methods that we learned about in this lecture -{doc}`Singular Value Decomposition ` to +## Overview + +This lecture applies computational methods that we learned about in the lecture +{doc}`Singular Value Decomposition ` to * first-order vector autoregressions (VARs) * dynamic mode decompositions (DMDs) -* connections between DMDs and first-order VARs +* connections between DMDs and first-order VARs -## First-Order Vector Autoregressions +{cite}`sargent2026dynamic` study these connections in detail. +We are especially interested in **tall and skinny** data sets in which the number of variables $m$ exceeds the number of time periods $n$. -We want to fit a **first-order vector autoregression** +Such data sets are common. -$$ -X_{t+1} = A X_t + C \epsilon_{t+1}, \quad \epsilon_{t+1} \perp X_t -$$ (eq:VARfirstorder) +For example, {cite}`SSY_CEX_2026` study quarterly Consumer Expenditure Survey data on 100 quantiles of each of three cross sections -- private income, post-tax-and-transfer income, and consumption -- together with aggregate income growth. -where $\epsilon_{t+1}$ is the time $t+1$ component of a sequence of i.i.d. $m \times 1$ random vectors with mean vector -zero and identity covariance matrix and where -the $ m \times 1 $ vector $ X_t $ is +That gives them $m = 301$ variables but only $n = 133$ quarterly observations. -$$ -X_t = \begin{bmatrix} X_{1,t} & X_{2,t} & \cdots & X_{m,t} \end{bmatrix}^\top -$$ (eq:Xvector) +They use a DMD to estimate a first-order VAR of rank $3$ for these 301 variables. -and where $\cdot ^\top $ again denotes complex transposition and $ X_{i,t} $ is variable $ i $ at time $ t $. +Along the way, we'll learn that +* a DMD computes a reduced-rank estimator of the coefficient matrix of a first-order VAR +* that estimator is a **principal components regression**: compress the cross section into a few principal components, then regress next period's data on them +* DMD **modes** are right eigenvectors of the estimated coefficient matrix +* matching **left** eigenvectors can be computed cheaply, and together with the modes they give an exact modal representation of the estimated VAR that is useful for forecasting and for computing long-run responses +We'll use the following imports. -We want to fit equation {eq}`eq:VARfirstorder`. +```{code-cell} ipython3 +import numpy as np +import matplotlib.pyplot as plt +``` +## First-order vector autoregressions -Our data are organized in an $ m \times (n+1) $ matrix $ \tilde X $ +We want to fit a **first-order vector autoregression** $$ -\tilde X = \begin{bmatrix} X_1 \mid X_2 \mid \cdots \mid X_n \mid X_{n+1} \end{bmatrix} +X_{t+1} = A X_t + C \epsilon_{t+1}, \quad \epsilon_{t+1} \perp X_t +$$ (eq:VARfirstorder) + +where $\epsilon_{t+1}$ is the time $t+1$ component of a sequence of i.i.d. $m \times 1$ random vectors with mean vector zero and identity covariance matrix and where the $m \times 1$ vector $X_t$ is + $$ +X_t = \begin{bmatrix} X_{1,t} & X_{2,t} & \cdots & X_{m,t} \end{bmatrix}^\top +$$ (eq:Xvector) + +Here $\cdot^\top$ denotes matrix transposition and $X_{i,t}$ is variable $i$ at time $t$. -where for $ t = 1, \ldots, n +1 $, the $ m \times 1 $ vector $ X_t $ is given by {eq}`eq:Xvector`. +Equation {eq}`eq:VARfirstorder` has no constant term. -Thus, we want to estimate a system {eq}`eq:VARfirstorder` that consists of $ m $ least squares regressions of **everything** on one lagged value of **everything**. +That is because we assume that each variable has already been transformed -- for example, by taking differences, by subtracting a common trend, and by subtracting its time-series mean -- so that $\{X_t\}$ can be regarded as a mean-zero, covariance-stationary process. -The $i$'th equation of {eq}`eq:VARfirstorder` is a regression of $X_{i,t+1}$ on the vector $X_t$. +For example, {cite}`SSY_CEX_2026` subtract a cross-section mean of log private income from every log quantile, use the growth rate of that mean as an additional variable, and then subtract time-series means from all variables. +Our data are organized in an $m \times (n+1)$ matrix $\tilde X$ -We proceed as follows. +$$ +\tilde X = \begin{bmatrix} X_1 \mid X_2 \mid \cdots \mid X_n \mid X_{n+1} \end{bmatrix} +$$ +where for $t = 1, \ldots, n+1$, the $m \times 1$ vector $X_t$ is given by {eq}`eq:Xvector`. -From $ \tilde X $, we form two $m \times n$ matrices +From $\tilde X$, we form two $m \times n$ matrices $$ X = \begin{bmatrix} X_1 \mid X_2 \mid \cdots \mid X_{n}\end{bmatrix} @@ -72,88 +90,79 @@ $$ X' = \begin{bmatrix} X_2 \mid X_3 \mid \cdots \mid X_{n+1}\end{bmatrix} $$ -Here $ ' $ is part of the name of the matrix $ X' $ and does not indicate matrix transposition. +Here $'$ is part of the name of the matrix $X'$ and does not indicate matrix transposition. -We use $\cdot^\top $ to denote matrix transposition or its extension to complex matrices. +In forming $X$ and $X'$, we have in each case dropped a column from $\tilde X$, the last column in the case of $X$, and the first column in the case of $X'$. -In forming $ X $ and $ X' $, we have in each case dropped a column from $ \tilde X $, the last column in the case of $ X $, and the first column in the case of $ X' $. +We want to estimate a system {eq}`eq:VARfirstorder` that consists of $m$ least squares regressions of **everything** on one lagged value of **everything**. -Evidently, $ X $ and $ X' $ are both $ m \times n $ matrices. +The $i$th equation of {eq}`eq:VARfirstorder` is a regression of $X_{i,t+1}$ on the vector $X_t$. -We denote the rank of $ X $ as $ p \leq \min(m, n) $. +So the $i$th row of an estimator $\hat A$ is a $1 \times m$ vector of regression coefficients of $X_{i,t+1}$ on $X_{j,t}, j = 1, \ldots, m$. -Two cases that interest us are +A least squares estimator $\hat A$ solves - * $ n > > m$, so that we have many more time series observations $n$ than variables $m$ - * $m > > n$, so that we have many more variables $m $ than time series observations $n$ +$$ +\min_{\check A} \| X' - \check A X \|_F +$$ (eq:ALSeqn) + +where $\| \cdot \|_F$ denotes the Frobenius norm of a matrix -At a general level that includes both of these special cases, a common formula describes the least squares estimator $\hat A$ of $A$. +$$ + \|B\|_F = \sqrt{ \sum_{i} \sum_{j} |B_{ij}|^2 } . +$$ -But important details differ. +We denote the rank of $X$ by $p \leq \min(m, n)$. -The common formula is +The **minimum-norm** solution of {eq}`eq:ALSeqn` is -$$ -\hat A = X' X^+ +$$ +\hat A = X' X^+ $$ (eq:commonA) -where $X^+$ is the pseudo-inverse of $X$. +where $X^+$ is the [Moore-Penrose pseudo-inverse](https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse) of $X$. -To read about the **Moore-Penrose pseudo-inverse** please see [Moore-Penrose pseudo-inverse](https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse) +Two cases interest us -Applicable formulas for the pseudo-inverse differ for our two cases. +* $n \gg m$, so that we have many more time series observations $n$ than variables $m$ +* $m \gg n$, so that we have many more variables $m$ than time series observations $n$ **Short-Fat Case:** -When $ n > > m$, so that we have many more time series observations $n$ than variables $m$ and when -$X$ has linearly independent **rows**, $X X^\top $ has an inverse and the pseudo-inverse $X^+$ is +When $n \gg m$ and $X$ has linearly independent **rows**, $X X^\top$ has an inverse and the pseudo-inverse $X^+$ is $$ -X^+ = X^\top (X X^\top )^{-1} +X^+ = X^\top (X X^\top )^{-1} $$ -Here $X^+$ is a **right-inverse** that verifies $ X X^+ = I_{m \times m}$. +Here $X^+$ is a **right-inverse** that verifies $X X^+ = I_{m \times m}$. -In this case, our formula {eq}`eq:commonA` for the least-squares estimator of the population matrix of regression coefficients $A$ becomes +In this case, problem {eq}`eq:ALSeqn` has a unique solution and formula {eq}`eq:commonA` becomes -$$ +$$ \hat A = X' X^\top (X X^\top )^{-1} $$ (eq:Ahatform101) +This formula for least-squares regression coefficients is widely used in econometrics to estimate vector autoregressions. -This formula for least-squares regression coefficients is widely used in econometrics. - -It is used to estimate vector autorgressions. - -The right side of formula {eq}`eq:Ahatform101` is proportional to the empirical cross second moment matrix of $X_{t+1}$ and $X_t$ times the inverse -of the second moment matrix of $X_t$. - - +The right side of formula {eq}`eq:Ahatform101` is proportional to the empirical cross second moment matrix of $X_{t+1}$ and $X_t$ times the inverse of the second moment matrix of $X_t$. **Tall-Skinny Case:** -When $m > > n$, so that we have many more attributes $m $ than time series observations $n$ and when $X$ has linearly independent **columns**, -$X^\top X$ has an inverse and the pseudo-inverse $X^+$ is +When $m \gg n$ and $X$ has linearly independent **columns**, $X^\top X$ has an inverse and the pseudo-inverse $X^+$ is $$ -X^+ = (X^\top X)^{-1} X^\top +X^+ = (X^\top X)^{-1} X^\top $$ -Here $X^+$ is a **left-inverse** that verifies $X^+ X = I_{n \times n}$. +Here $X^+$ is a **left-inverse** that verifies $X^+ X = I_{n \times n}$. -In this case, our formula {eq}`eq:commonA` for a least-squares estimator of $A$ becomes +In this case, formula {eq}`eq:commonA` becomes $$ -\hat A = X' (X^\top X)^{-1} X^\top +\hat A = X' (X^\top X)^{-1} X^\top $$ (eq:hatAversion0) -Please compare formulas {eq}`eq:Ahatform101` and {eq}`eq:hatAversion0` for $\hat A$. - -Here we are especially interested in formula {eq}`eq:hatAversion0`. - -The $ i $th row of $ \hat A $ is an $ m \times 1 $ vector of regression coefficients of $ X_{i,t+1} $ on $ X_{j,t}, j = 1, \ldots, m $. - - If we use formula {eq}`eq:hatAversion0` to calculate $\hat A X$ we find that $$ @@ -164,677 +173,804 @@ so that the regression equation **fits perfectly**. This is a typical outcome in an **underdetermined least-squares** model. +In fact, every matrix $\check A = \hat A + N$ with $N X = 0$ also fits perfectly, so problem {eq}`eq:ALSeqn` has infinitely many solutions. -To reiterate, in the **tall-skinny** case (described in {doc}`Singular Value Decomposition `) in which we have a number $n$ of observations that is small relative to the number $m$ of -attributes that appear in the vector $X_t$, we want to fit equation {eq}`eq:VARfirstorder`. +Formula {eq}`eq:hatAversion0` selects the one with the smallest Frobenius norm. -We confront the facts that the least squares estimator is underdetermined and that the regression equation fits perfectly. +A perfect in-sample fit is a warning sign: with more coefficients than observations, the estimator $\hat A$ is tracking noise as well as signal. +We'll respond by constructing a **reduced-rank** estimator that no longer fits perfectly. -To proceed, we'll want efficiently to calculate the pseudo-inverse $X^+$. +### Computing $X^+$ with an SVD -The pseudo-inverse $X^+$ will be a component of our estimator of $A$. +An efficient way to compute the pseudo-inverse $X^+$ is to start with a **reduced** singular value decomposition -As our estimator $\hat A$ of $A$ we want to form an $m \times m$ matrix that solves the least-squares best-fit problem +$$ +X = U \Sigma V^\top +$$ (eq:SVDDMD) -$$ -\hat A = \textrm{argmin}_{\check A} || X' - \check A X ||_F -$$ (eq:ALSeqn) +where $U$ is an $m \times p$ matrix, $\Sigma$ is a $p \times p$ diagonal matrix of positive singular values $\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_p > 0$, and $V$ is an $n \times p$ matrix. -where $|| \cdot ||_F$ denotes the Frobenius (or Euclidean) norm of a matrix. +As we saw in {doc}`Singular Value Decomposition `, for a reduced SVD -The Frobenius norm is defined as +* $U^\top U = I_{p \times p}$ and $V^\top V = I_{p \times p}$, but +* $U U^\top \neq I_{m \times m}$ unless $p = m$, and $V V^\top \neq I_{n \times n}$ unless $p = n$. -$$ - ||A||_F = \sqrt{ \sum_{i=1}^m \sum_{j=1}^m |A_{ij}|^2 } -$$ +For any rank $p$, the Moore-Penrose pseudo-inverse of $X$ is +$$ +X^{+} = V \Sigma^{-1} U^\top +$$ (eq:Xplusformula) -The minimizer of the right side of equation {eq}`eq:ALSeqn` is +When $X$ has linearly independent columns, so that $p = n$ and $V$ is a square orthogonal matrix, we can confirm that formula {eq}`eq:Xplusformula` agrees with our earlier formula: $$ -\hat A = X' X^{+} -$$ (eq:hatAform) - -where the (possibly huge) $ n \times m $ matrix $ X^{+} = (X^\top X)^{-1} X^\top $ is again a pseudo-inverse of $ X $. +\begin{aligned} +(X^\top X)^{-1} X^\top + & = (V \Sigma U^\top U \Sigma V^\top )^{-1} V \Sigma U^\top \\ + & = (V \Sigma^2 V^\top )^{-1} V \Sigma U^\top \\ + & = V \Sigma^{-2} V^\top V \Sigma U^\top \\ + & = V \Sigma^{-1} U^\top +\end{aligned} +$$ +where the third line uses $V^{-1} = V^\top$. +Substituting formula {eq}`eq:Xplusformula` into formula {eq}`eq:commonA` gives +$$ +\hat A = X' V \Sigma^{-1} U^\top +$$ (eq:AhatSVDformula) -For some situations that we are interested in, $X^\top X $ can be close to singular, a situation that makes some numerical algorithms be inaccurate. +Formula {eq}`eq:AhatSVDformula` shows a weakness of $\hat A$. -To acknowledge that possibility, we'll use efficient algorithms to constructing -a **reduced-rank approximation** of $\hat A$ in formula {eq}`eq:hatAversion0`. +Small singular values $\sigma_j$ enter through $1/\sigma_j$. -Such an approximation to our vector autoregression will no longer fit perfectly. - +So directions in which the data $X$ vary very little receive very large weights, which amplifies noise. -The $ i $th row of $ \hat A $ is an $ m \times 1 $ vector of regression coefficients of $ X_{i,t+1} $ on $ X_{j,t}, j = 1, \ldots, m $. +## A reduced-rank VAR -An efficient way to compute the pseudo-inverse $X^+$ is to start with a singular value decomposition +A remedy is to keep only the $r < p$ largest singular values. +Let $U_r$ be the $m \times r$ matrix consisting of the first $r$ columns of $U$, let $\Sigma_r$ be the $r \times r$ upper-left block of $\Sigma$, and let $V_r$ be the $n \times r$ matrix consisting of the first $r$ columns of $V$. +Our **rank-$r$ estimator** of $A$ is $$ -X = U \Sigma V^\top -$$ (eq:SVDDMD) +\hat A_r = X' V_r \Sigma_r^{-1} U_r^\top +$$ (eq:Ahat_r) -where we remind ourselves that for a **reduced** SVD, $X$ is an $m \times n$ matrix of data, $U$ is an $m \times p$ matrix, $\Sigma$ is a $p \times p$ matrix, and $V$ is an $n \times p$ matrix. +The $m \times m$ matrix $\hat A_r$ has rank at most $r$. -We can efficiently construct the pertinent pseudo-inverse $X^+$ -by recognizing the following string of equalities. +Notice that $V_r \Sigma_r^{-1} U_r^\top$ is the pseudo-inverse of $\hat X_r = U_r \Sigma_r V_r^\top$, the best rank-$r$ approximation of $X$ described by the Eckart-Young theorem in {doc}`Singular Value Decomposition `. -$$ -\begin{aligned} -X^{+} & = (X^\top X)^{-1} X^\top \\ - & = (V \Sigma U^\top U \Sigma V^\top )^{-1} V \Sigma U^\top \\ - & = (V \Sigma \Sigma V^\top )^{-1} V \Sigma U^\top \\ - & = V \Sigma^{-1} \Sigma^{-1} V^\top V \Sigma U^\top \\ - & = V \Sigma^{-1} U^\top -\end{aligned} -$$ (eq:efficientpseudoinverse) +So $\hat A_r = X' \hat X_r^+$: we obtain the rank-$r$ estimator by replacing $X$ with its Eckart-Young approximation in the least squares formula $\hat A = X' X^+$. +This is the estimator that a **dynamic mode decomposition** computes. -(Since we are in the $m > > n$ case in which $V^\top V = I_{p \times p}$ in a reduced SVD, we can use the preceding -string of equalities for a reduced SVD as well as for a full SVD.) +Dynamic mode decomposition was introduced by {cite}`schmid2010`. -Thus, we shall construct a pseudo-inverse $ X^+ $ of $ X $ by using -a singular value decomposition of $X$ in equation {eq}`eq:SVDDMD` to compute +You can read about it in {cite}`DMD_book` and {cite}`Brunton_Kutz_2019` (section 7.2). +Associated with $\hat A_r$ are fitted residuals and their covariance matrix $$ -X^{+} = V \Sigma^{-1} U^\top -$$ (eq:Xplusformula) +\hat a_{t+1} = X_{t+1} - \hat A_r X_t, \qquad +\hat \Omega = \frac{1}{n} \sum_{t=1}^{n} \hat a_{t+1} \hat a_{t+1}^\top += \frac{1}{n} (X' - \hat A_r X)(X' - \hat A_r X)^\top +$$ (eq:residualsOmega) -where the matrix $ \Sigma^{-1} $ is constructed by replacing each non-zero element of $ \Sigma $ with $ \sigma_j^{-1} $. +Thus, our estimated rank-$r$ VAR is -We can use formula {eq}`eq:Xplusformula` together with formula {eq}`eq:hatAform` to compute the matrix $ \hat A $ of regression coefficients. +$$ +X_{t+1} = \hat A_r X_t + \hat a_{t+1} +$$ (eq:reducedVAR) -Thus, our estimator $\hat A = X' X^+$ of the $m \times m$ matrix of coefficients $A$ is +With $r < p$, the regression no longer fits perfectly and $\hat \Omega$ is not zero. -$$ -\hat A = X' V \Sigma^{-1} U^\top -$$ (eq:AhatSVDformula) +In applications, $r$ is often very small, e.g., three or fewer. +A common way to choose $r$ is to plot the singular values $\sigma_1, \sigma_2, \ldots$ and to truncate where they stop falling rapidly and level off -- an "elbow rule." +{cite}`SSY_CEX_2026` used this rule to choose $r = 3$. -## Dynamic Mode Decomposition (DMD) +### Interpretation as a principal components regression +We can interpret formula {eq}`eq:Ahat_r` statistically. +Define the $r \times 1$ vector -We turn to the $ m >>n $ **tall and skinny** case associated with **Dynamic Mode Decomposition**. +$$ +f_t = U_r^\top X_t +$$ -Here an $ m \times n+1 $ data matrix $ \tilde X $ contains many more attributes (or variables) $ m $ than time periods $ n+1 $. +whose components are the first $r$ **principal components** of the cross section $X_t$ (see the discussion of principal components analysis in {doc}`Singular Value Decomposition `). +Stacking them gives the $r \times n$ matrix $F = U_r^\top X = \Sigma_r V_r^\top$. -Dynamic mode decomposition was introduced by {cite}`schmid2010`, +The least-squares regression of $X_{t+1}$ on $f_t$ has coefficient matrix + +$$ +X' F^\top (F F^\top)^{-1} = X' V_r \Sigma_r (\Sigma_r^2)^{-1} = X' V_r \Sigma_r^{-1} +$$ -You can read about Dynamic Mode Decomposition {cite}`DMD_book` and {cite}`Brunton_Kutz_2019` (section 7.2). +so that +$$ +\hat A_r X_t = \left( X' V_r \Sigma_r^{-1} \right) f_t . +$$ -**Dynamic Mode Decomposition** (DMD) computes a rank $ r < p $ approximation to the least squares regression coefficients $ \hat A $ described by formula {eq}`eq:AhatSVDformula`. +Similarly, the least-squares regression of $f_{t+1} = U_r^\top X_{t+1}$ on $f_t$ has $r \times r$ coefficient matrix - -We'll build up gradually to a formulation that is useful in applications. +$$ +\tilde A = U_r^\top X' V_r \Sigma_r^{-1} +$$ (eq:Atilde) +So $\tilde A$ is the first-order VAR coefficient matrix for the $r$ principal components. -We'll do this by describing three alternative representations of our first-order linear dynamic system, i.e., our vector autoregression. +Note that $\tilde A = U_r^\top \hat A_r U_r$. -**Guide to three representations:** In practice, we'll mainly be interested in Representation 3. +Thus, a DMD proceeds in two steps: -We use the first two representations to present some useful intermediate steps that help us to appreciate what is under the hood of Representation 3. +1. it compresses the cross section $X_t$ into $r$ principal components $f_t$, using only the variation in $X$ +1. it regresses next period's data $X_{t+1}$ on those principal components -In applications, we'll use only a small subset of **DMD modes** to approximate dynamics. +Statisticians call this a **principal components regression**. -We use such a small subset of DMD modes to construct a reduced-rank approximation to $A$. +It resembles the "diffusion index" forecasting method of {cite}`stock_watson2002`. -To do that, we'll want to use the **reduced** SVD's affiliated with representation 3, not the **full** SVD's affiliated with representations 1 and 2. +```{note} +The estimator $\hat A_r$ is **not** the matrix of rank $r$ that minimizes $\| X' - \check A X \|_F$. +That minimizer is the **reduced-rank regression** estimator of {cite}`anderson1951`. -**Guide to impatient reader:** In our applications, we'll be using Representation 3. +It chooses the $r$ directions that best *predict* $X_{t+1}$, while $\hat A_r$ uses the $r$ directions that best *describe* $X_t$. -You might want to skip the stage-setting representations 1 and 2 on first reading. +{ref}`var_dmd_ex2` compares the two estimators. +``` -+++ +## The DMD algorithm -## Representation 1 - -In this representation, we shall use a **full** SVD of $X$. +Here is the DMD algorithm. -We use the $m$ **columns** of $U$, and thus the $m$ **rows** of $U^\top $, to define a $m \times 1$ vector $\tilde b_t$ as +1. Compute the reduced SVD $X = U \Sigma V^\top$. +1. Retain the $r$ largest singular values to form $U_r$, $\Sigma_r$, and $V_r$. +1. Form the $r \times r$ matrix $\tilde A = U_r^\top X' V_r \Sigma_r^{-1}$ and compute its eigendecomposition + $$ + \tilde A W = W \Lambda + $$ (eq:tildeAeigenred) -$$ -\tilde b_t = U^\top X_t . -$$ (eq:tildeXdef2) + where $\Lambda$ is a diagonal matrix of eigenvalues $\lambda_1, \ldots, \lambda_r$ and the columns of $W$ are corresponding eigenvectors. +1. Form the $m \times r$ matrix $\Phi$ and the $r \times m$ matrix $\Psi$ -The original data $X_t$ can be represented as + $$ + \Phi = X' V_r \Sigma_r^{-1} W, \qquad \Psi = (W \Lambda)^{-1} U_r^\top + $$ (eq:PhiPsi) -$$ -X_t = U \tilde b_t -$$ (eq:Xdecoder) +The eigenvalues $\lambda_i$ on the diagonal of $\Lambda$ are **DMD eigenvalues** and the columns $\phi_i$ of $\Phi$ are **DMD modes**. -(Here we use $b$ to remind ourselves that we are creating a **basis** vector.) +We'll see below that the rows $\psi_i$ of $\Psi$ are matching **left** eigenvectors. -Since we are now using a **full** SVD, $U U^\top = I_{m \times m}$. +Notice that the heavy lifting consists of an SVD of the $m \times n$ matrix $X$ and an eigendecomposition of the small $r \times r$ matrix $\tilde A$. -So it follows from equation {eq}`eq:tildeXdef2` that we can reconstruct $X_t$ from $\tilde b_t$. +We never have to form or decompose an $m \times m$ matrix. -In particular, +Throughout, we assume that $\tilde A$ has $r$ linearly independent eigenvectors and nonzero eigenvalues, which is the generic case. +The matrix $\tilde A$ is real but not symmetric, so some eigenvalues can come in complex-conjugate pairs. +In that case, the corresponding columns of $\Phi$ and rows of $\Psi$ also come in complex-conjugate pairs, while products such as $\Phi \Lambda^j \Psi$ are real. - * Equation {eq}`eq:tildeXdef2` serves as an **encoder** that **rotates** the $m \times 1$ vector $X_t$ to become an $m \times 1$ vector $\tilde b_t$ - - * Equation {eq}`eq:Xdecoder` serves as a **decoder** that **reconstructs** the $m \times 1$ vector $X_t$ by rotating the $m \times 1$ vector $\tilde b_t$ +For complex matrices, we use $\cdot^*$ to denote the conjugate transpose. +```{note} +{cite}`tu_Rowley` scale the modes differently, defining them as the columns of $\Phi \Lambda^{-1}$. +With that scaling, the matching left eigenvectors are the rows of $\Lambda \Psi = W^{-1} U_r^\top$. -Define a transition matrix for an $m \times 1$ basis vector $\tilde b_t$ by +Rescaling a mode and its left eigenvector in opposite directions changes none of the formulas below. +``` -$$ -\tilde A = U^\top \hat A U -$$ (eq:Atilde0) +## Right eigenvectors: DMD modes -We can recover $\hat A$ from +The following result of {cite}`tu_Rowley` explains why the columns of $\Phi$ are called modes. -$$ -\hat A = U \tilde A U^\top -$$ +```{prf:proposition} +:label: prop-dmd-eigenvectors -Dynamics of the $m \times 1$ basis vector $\tilde b_t$ are governed by +The columns of $\Phi$ are eigenvectors of $\hat A_r$: $$ -\tilde b_{t+1} = \tilde A \tilde b_t +\hat A_r \Phi = \Phi \Lambda $$ +``` -To construct forecasts $\overline X_t$ of future values of $X_t$ conditional on $X_1$, we can apply decoders (i.e., rotators) to both sides of this equation and deduce +```{prf:proof} +Using formulas {eq}`eq:Ahat_r`, {eq}`eq:Atilde`, and {eq}`eq:PhiPsi`, $$ -\overline X_{t+1} = U \tilde A^t U^\top X_1 +\begin{aligned} + \hat A_r \Phi & = (X' V_r \Sigma_r^{-1} U_r^\top) (X' V_r \Sigma_r^{-1} W) \\ + & = X' V_r \Sigma_r^{-1} \tilde A W \\ + & = X' V_r \Sigma_r^{-1} W \Lambda \\ + & = \Phi \Lambda +\end{aligned} $$ +``` -where we use $\overline X_{t+1}, t \geq 1 $ to denote a forecast. - -+++ - -## Representation 2 - +Equating columns gives $\hat A_r \phi_i = \lambda_i \phi_i$ for $i = 1, \ldots, r$. -This representation is related to one originally proposed by {cite}`schmid2010`. +We can say more. -It can be regarded as an intermediate step on the way to obtaining a related representation 3 to be presented later +Write $\hat A_r = B U_r^\top$, where $B = X' V_r \Sigma_r^{-1}$ is $m \times r$. +Because the nonzero eigenvalues of $B U_r^\top$ coincide with those of $U_r^\top B = \tilde A$, the eigenvalues of $\hat A_r$ are $\lambda_1, \ldots, \lambda_r$ together with $m - r$ zeros. -As with Representation 1, we continue to +So the estimated VAR {eq}`eq:reducedVAR` is stable if and only if $\max_i |\lambda_i| < 1$, a condition we can check with an $r \times r$ matrix. -* use a **full** SVD and **not** a reduced SVD +(Also see {cite}`DDSE_book`, p. 238.) +### Projected modes +{cite}`schmid2010` originally worked with the **projected modes** $U_r W$ instead of $\Phi$. -As we observed and illustrated in a lecture about the {doc}`Singular Value Decomposition ` +Because $U_r^\top \Phi = \tilde A W = W \Lambda$, we have - * (a) for a full SVD $U U^\top = I_{m \times m} $ and $U^\top U = I_{p \times p}$ are both identity matrices - - * (b) for a reduced SVD of $X$, $U^\top U $ is not an identity matrix. +$$ +U_r U_r^\top \Phi = U_r W \Lambda +$$ -As we shall see later, a full SVD is too confining for what we ultimately want to do, namely, cope with situations in which $U^\top U$ is **not** an identity matrix because we use a reduced SVD of $X$. +so, up to scale, the projected modes are orthogonal projections of the modes $\Phi$ onto the column space of $U_r$. -But for now, let's proceed under the assumption that we are using a full SVD so that requirements (a) and (b) are both satisfied. +The projected modes are eigenvectors of $\hat A_r$ only when $\Phi = U_r W \Lambda$, that is, when the columns of $\Phi$ lie in the column space of $U_r$. - +That need not be true, because $\hat A_r U_r W = \Phi$. -Form an eigendecomposition of the $m \times m$ matrix $\tilde A = U^\top \hat A U$ defined in equation {eq}`eq:Atilde0`: +The construction $\Phi = \hat A_r U_r W$ of {cite}`tu_Rowley` delivers eigenvectors of $\hat A_r$ in all cases. -$$ -\tilde A = W \Lambda W^{-1} -$$ (eq:tildeAeigen) +## Left eigenvectors and a modal representation -where $\Lambda$ is a diagonal matrix of eigenvalues and $W$ is an $m \times m$ -matrix whose columns are eigenvectors corresponding to rows (eigenvalues) in -$\Lambda$. +Now we come to the rows of $\Psi$. -When $U U^\top = I_{m \times m}$, as is true with a full SVD of $X$, it follows that +```{prf:proposition} +:label: prop-dmd-left -$$ -\hat A = U \tilde A U^\top = U W \Lambda W^{-1} U^\top -$$ (eq:eqeigAhat) +Let $\Psi = (W \Lambda)^{-1} U_r^\top$. Then -According to equation {eq}`eq:eqeigAhat`, the diagonal matrix $\Lambda$ contains eigenvalues of $\hat A$ and corresponding eigenvectors of $\hat A$ are columns of the matrix $UW$. +1. $\Psi \Phi = I_{r \times r}$ +1. $\Psi \hat A_r = \Lambda \Psi$, so that the rows of $\Psi$ are left eigenvectors of $\hat A_r$ +1. $\hat A_r = \Phi \Lambda \Psi = \sum_{i=1}^r \lambda_i \phi_i \psi_i$ +``` -It follows that the systematic (i.e., not random) parts of the $X_t$ dynamics captured by our first-order vector autoregressions are described by +```{prf:proof} +For the first claim, $$ -X_{t+1} = U W \Lambda W^{-1} U^\top X_t +\Psi \Phi = \Lambda^{-1} W^{-1} U_r^\top X' V_r \Sigma_r^{-1} W += \Lambda^{-1} W^{-1} \tilde A W = \Lambda^{-1} \Lambda = I . $$ -Multiplying both sides of the above equation by $W^{-1} U^\top $ gives +For the second claim, since $U_r^\top \hat A_r = U_r^\top X' V_r \Sigma_r^{-1} U_r^\top = \tilde A U_r^\top$, -$$ -W^{-1} U^\top X_{t+1} = \Lambda W^{-1} U^\top X_t $$ - -or - -$$ -\hat b_{t+1} = \Lambda \hat b_t -$$ - -where our **encoder** is - -$$ -\hat b_t = W^{-1} U^\top X_t +\Psi \hat A_r = \Lambda^{-1} W^{-1} \tilde A U_r^\top += \Lambda^{-1} \Lambda W^{-1} U_r^\top = W^{-1} U_r^\top = \Lambda \Psi . $$ -and our **decoder** is +For the third claim, $$ -X_t = U W \hat b_t +\Phi \Lambda \Psi = X' V_r \Sigma_r^{-1} W \Lambda \Lambda^{-1} W^{-1} U_r^\top += X' V_r \Sigma_r^{-1} U_r^\top = \hat A_r . $$ +``` -We can use this representation to construct a predictor $\overline X_{t+1}$ of $X_{t+1}$ conditional on $X_1$ via: +Define the $r \times 1$ vector of **modal coordinates** $$ -\overline X_{t+1} = U W \Lambda^t W^{-1} U^\top X_1 -$$ (eq:DSSEbookrepr) - - -In effect, -{cite}`schmid2010` defined an $m \times m$ matrix $\Phi_s$ as +b_t = \Psi X_t = (W \Lambda)^{-1} U_r^\top X_t +$$ (eq:modalcoords) -$$ -\Phi_s = UW -$$ (eq:Phisfull) +Computing $b_t$ requires only the principal components $f_t = U_r^\top X_t$ and an $r \times r$ matrix. -and a generalized inverse +Because $\Psi \Phi = I$, iterating $\hat A_r = \Phi \Lambda \Psi$ gives $$ -\Phi_s^+ = W^{-1}U^\top -$$ (eq:Phisfullinv) +\hat A_r^j = \Phi \Lambda^j \Psi, \quad j \geq 1 +$$ (eq:Ahatpower) -{cite}`schmid2010` then represented equation {eq}`eq:DSSEbookrepr` as +Consequently, the estimated VAR {eq}`eq:reducedVAR` implies $$ -\overline X_{t+1} = \Phi_s \Lambda^t \Phi_s^+ X_1 -$$ (eq:schmidrep) - -Components of the basis vector $ \hat b_t = W^{-1} U^\top X_t \equiv \Phi_s^+ X_t$ are -DMD **projected modes**. - -To understand why they are called **projected modes**, notice that - -$$ -\Phi_s^+ = ( \Phi_s^\top \Phi_s)^{-1} \Phi_s^\top +b_{t+1} = \Lambda b_t + \Psi \hat a_{t+1} $$ -so that the $m \times p$ matrix +and the $j$-step-ahead forecast of $X_{t+j}$ conditional on $X_t$ is $$ -\hat b = \Phi_s^+ X -$$ - -is a matrix of regression coefficients of the $m \times n$ matrix $X$ on the $m \times p$ matrix $\Phi_s$. - -We'll say more about this interpretation in a related context when we discuss representation 3, which was suggested by Tu et al. {cite}`tu_Rowley`. +\hat A_r^j X_t = \Phi \Lambda^j b_t = \sum_{i=1}^r \phi_i \lambda_i^j b_{i,t} +$$ (eq:modalforecast) -It is more appropriate to use representation 3 when, as is often the case in practice, we want to use a reduced SVD. +Each mode $\phi_i$ contributes a component whose amplitude decays geometrically at rate $\lambda_i$. +The matrix $\Phi \Psi$ is an oblique projection onto the column space of $\Phi$. +It splits $X_t$ into $\Phi b_t$ plus a remainder $(I - \Phi \Psi) X_t$ that $\hat A_r$ sends to zero, because $\hat A_r (I - \Phi \Psi) = \Phi \Lambda \Psi - \Phi \Lambda \Psi \Phi \Psi = 0$. +### A tempting alternative that does not work -## Representation 3 +The Moore-Penrose pseudo-inverse $\Phi^+ = (\Phi^* \Phi)^{-1} \Phi^*$ also satisfies $\Phi^+ \Phi = I$. -Departing from the procedures used to construct Representations 1 and 2, each of which deployed a **full** SVD, we now use a **reduced** SVD. +The vector $\check b_t = \Phi^+ X_t$ is a vector of least squares regression coefficients of $X_t$ on the columns of $\Phi$, and $\Phi \check b_t$ is the orthogonal projection of $X_t$ onto the column space of $\Phi$. -Again, we let $p \leq \textrm{min}(m,n)$ be the rank of $X$. +That makes $\check b_t$ a natural **descriptive** summary of how much each mode is present in $X_t$. -Construct a **reduced** SVD +But in general $$ -X = \tilde U \tilde \Sigma \tilde V^\top , +\Phi \Lambda \Phi^+ \neq \hat A_r $$ -where now $\tilde U$ is $m \times p$, $\tilde \Sigma$ is $ p \times p$, and $\tilde V^\top $ is $p \times n$. +so $\check b_t$ does **not** evolve according to $\check b_{t+1} = \Lambda \check b_t$ under the estimated VAR, and $\Phi \Lambda^j \Phi^+ X_t$ is not the VAR forecast $\hat A_r^j X_t$. -Our minimum-norm least-squares approximator of $A$ now has representation +The reason is that $\Phi^+$ and $\Psi$ are both left inverses of $\Phi$, but their rows span different spaces: -$$ -\hat A = X' \tilde V \tilde \Sigma^{-1} \tilde U^\top -$$ (eq:Ahatwithtildes) +* the rows of $\Phi^+$ span the same space as the rows of $\Phi^*$ +* the rows of $\Psi$ span the same space as the rows of $U_r^\top$ +The two left inverses coincide if and only if the column spaces of $\Phi$ and $U_r$ are the same. -**Computing Dominant Eigenvectors of $\hat A$** +The column space of $\Phi$ is the column space of $X' V_r$, which describes where the data go **next** period. -We begin by paralleling a step used to construct Representation 1, define a transition matrix for a rotated $p \times 1$ state $\tilde b_t$ by +The column space of $U_r$ is spanned by the $r$ directions in which the **current** data $X$ vary most. -$$ -\tilde A =\tilde U^\top \hat A \tilde U -$$ (eq:Atildered) +These two spaces coincide in special cases, for example when noise-free data are driven by exactly $r$ factors (see {ref}`var_dmd_ex1`), but not in general. +```{note} +The formula $(W \Lambda)^{-1} U_r^\top X_1$ appears in {cite}`DDSE_book` (p. 240), where it is presented as a computationally cheap approximation to $\Phi^+ X_1$. -**Interpretation as projection coefficients** +{prf:ref}`prop-dmd-left` shows that it is exactly the vector of coordinates that is consistent with the dynamics of the estimated VAR. +``` +### A state-space representation -{cite}`DDSE_book` remark that $\tilde A$ can be interpreted in terms of a projection of $\hat A$ onto the $p$ modes in $\tilde U$. +The modal representation lets us write the estimated VAR as a linear state-space system with an $r \times 1$ state vector. -To verify this, first note that, because $ \tilde U^\top \tilde U = I$, it follows that +Define $$ -\tilde A = \tilde U^\top \hat A \tilde U = \tilde U^\top X' \tilde V \tilde \Sigma^{-1} \tilde U^\top \tilde U -= \tilde U^\top X' \tilde V \tilde \Sigma^{-1} \tilde U^\top -$$ (eq:tildeAverify) +\hat x_t = \Lambda \Psi X_{t-1} = W^{-1} U_r^\top X_{t-1} = W^{-1} f_{t-1} +$$ +Then $X_t = \hat A_r X_{t-1} + \hat a_t = \Phi \Lambda \Psi X_{t-1} + \hat a_t$ and $\hat x_{t+1} = \Lambda \Psi X_t$ imply - +$$ +\begin{aligned} +\hat x_{t+1} & = \Lambda \hat x_t + \Lambda \Psi \hat a_t \\ +X_t & = \Phi \hat x_t + \hat a_t +\end{aligned} +$$ (eq:dmd_statespace) -Next, we'll just compute the regression coefficients in a projection of $\hat A$ on $\tilde U$ using a standard least-squares formula +where we have used $\Psi \Phi = I$. -$$ -(\tilde U^\top \tilde U)^{-1} \tilde U^\top \hat A = (\tilde U^\top \tilde U)^{-1} \tilde U^\top X' \tilde V \tilde \Sigma^{-1} \tilde U^\top = -\tilde U^\top X' \tilde V \tilde \Sigma^{-1} \tilde U^\top = \tilde A . -$$ +In representation {eq}`eq:dmd_statespace` -Thus, we have verified that $\tilde A$ is a least-squares projection of $\hat A$ onto $\tilde U$. +* the state $\hat x_t$ is the vector of lagged principal components expressed in the eigenvector basis of the principal-component VAR $\tilde A$ +* the transition matrix $\Lambda$ is diagonal, so the $r$ components of the state evolve independently, apart from correlated shocks +* the columns of $\Phi$ are **loadings** of the $m$ observed variables on the $r$ components of the state +* a single $m \times 1$ vector of residuals $\hat a_t$ drives both equations; by the least squares normal equations, it is orthogonal in sample to $f_{t-1}$ and therefore to $\hat x_t$ -**An Inverse Challenge** +{cite}`sargent2026dynamic` connect DMDs to linear state-space models of this form. +## Long-run responses -Because we are using a reduced SVD, $\tilde U \tilde U^\top \neq I$. +Modal representation {eq}`eq:Ahatpower` also simplifies calculations of long-run responses. -Consequently, +If $|\lambda_i| < 1$ for all $i$, then $$ -\hat A \neq \tilde U \tilde A \tilde U^\top , +\sum_{j=1}^\infty \hat A_r^j = \Phi \Lambda (I - \Lambda)^{-1} \Psi $$ -so we can't simply recover $\hat A$ from $\tilde A$ and $\tilde U$. - -**A Blind Alley** - -We can start by hoping for the best and proceeding to construct an eigendecomposition of the $p \times p$ matrix $\tilde A$: +and $$ - \tilde A = \tilde W \Lambda \tilde W^{-1} -$$ (eq:tildeAeigenred) +(I - \hat A_r)^{-1} = I + \Phi \Lambda (I - \Lambda)^{-1} \Psi +$$ (eq:longrun) -where $\Lambda$ is a diagonal matrix of $p$ eigenvalues and the columns of $\tilde W$ -are corresponding eigenvectors. +The matrix $\Lambda (I - \Lambda)^{-1}$ is diagonal with entries $\lambda_i / (1 - \lambda_i)$. +So long-run responses are dominated by the most persistent modes. -Mimicking our procedure in Representation 2, we cross our fingers and compute an $m \times p$ matrix +For example, the rank-3 VAR of {cite}`SSY_CEX_2026` has DMD eigenvalues $0.963$, $0.865$, and $0.547$, with long-run weights of roughly $26$, $6.4$, and $1.2$. -$$ -\tilde \Phi_s = \tilde U \tilde W -$$ (eq:Phisred) +Formulas like {eq}`eq:longrun` appear in **additive functionals**, which are studied in {doc}`advanced:additive_functionals`. -that corresponds to {eq}`eq:Phisfull` for a full SVD. +Suppose that the first component of $X_t$ is the growth rate of the logarithm $Y_t$ of some aggregate, measured as a deviation from its mean $\nu$, so that $Y_{t+1} - Y_t - \nu = e_1 X_{t+1}$ where $e_1 = \begin{bmatrix} 1 & 0 & \cdots & 0 \end{bmatrix}$. -At this point, where $\hat A$ is given by formula {eq}`eq:Ahatwithtildes` it is interesting to compute $\hat A \tilde \Phi_s$: +If $X_t$ obeys the VAR {eq}`eq:reducedVAR`, then $Y_t$ can be decomposed into a deterministic trend $t \nu$, a martingale $\sum_{j=1}^t H \hat a_j$, and a stationary component $-g X_t$, plus a constant, where $$ -\begin{aligned} -\hat A \tilde \Phi_s & = (X' \tilde V \tilde \Sigma^{-1} \tilde U^\top ) (\tilde U \tilde W) \\ - & = X' \tilde V \tilde \Sigma^{-1} \tilde W \\ - & \neq (\tilde U \tilde W) \Lambda \\ - & = \tilde \Phi_s \Lambda - \end{aligned} +g = e_1 \hat A_r (I - \hat A_r)^{-1}, \qquad H = e_1 (I - \hat A_r)^{-1} = e_1 + g $$ - -That -$ \hat A \tilde \Phi_s \neq \tilde \Phi_s \Lambda $ means that, unlike the corresponding situation in Representation 2, columns of $\tilde \Phi_s = \tilde U \tilde W$ -are **not** eigenvectors of $\hat A$ corresponding to eigenvalues on the diagonal of matix $\Lambda$. -**An Approach That Works** +{cite}`SSY_CEX_2026` use this decomposition for aggregate income and for all 300 CEX quantiles. -Continuing our quest for eigenvectors of $\hat A$ that we **can** compute with a reduced SVD, let's define an $m \times p$ matrix -$\Phi$ as +Formula {eq}`eq:longrun` implies $$ -\Phi \equiv \hat A \tilde \Phi_s = X' \tilde V \tilde \Sigma^{-1} \tilde W -$$ (eq:Phiformula) - -It turns out that columns of $\Phi$ **are** eigenvectors of $\hat A$. - -This is a consequence of a result established by Tu et al. {cite}`tu_Rowley` that we now present. - - - - - -```{prf:proposition} -:label: prop-dmd-eigenvectors +g = e_1 \Phi \Lambda (I - \Lambda)^{-1} \Psi +$$ -The $p$ columns of $\Phi$ are eigenvectors of $\hat A$. -``` +so computing $g$ and $H$ requires no inversion of an $m \times m$ matrix. -```{prf:proof} -From formula {eq}`eq:Phiformula` we have +## A simulated example -$$ -\begin{aligned} - \hat A \Phi & = (X' \tilde V \tilde \Sigma^{-1} \tilde U^\top ) (X' \tilde V \Sigma^{-1} \tilde W) \cr - & = X' \tilde V \tilde \Sigma^{-1} \tilde A \tilde W \cr - & = X' \tilde V \tilde \Sigma^{-1}\tilde W \Lambda \cr - & = \Phi \Lambda - \end{aligned} -$$ +To see these objects at work, we simulate a tall-skinny data set that shares some features with the CEX data of {cite}`SSY_CEX_2026`. -so that +There are $m = 300$ cross-section variables, which we can think of as quantiles indexed by $q \in [0, 1]$, and $n + 1 = 134$ quarterly observations. -$$ -\hat A \Phi = \Phi \Lambda . -$$ (eq:APhiLambda) +The cross section is driven by three latent factors $z_t$ that obey a first-order VAR with a diagonal transition matrix -Let $\phi_i$ be the $i$th column of $\Phi$ and $\lambda_i$ be the corresponding $i$ eigenvalue of $\tilde A$ from decomposition {eq}`eq:tildeAeigenred`. +$$ +z_{t+1} = D z_t + w_{t+1}, \qquad D = \textrm{diag}(0.95, 0.80, 0.50) +$$ -Equating the $m \times 1$ vectors that appear on the two sides of equation {eq}`eq:APhiLambda` gives +and the observed data are $$ -\hat A \phi_i = \lambda_i \phi_i . +X_t = G z_t + \sigma_v v_t $$ -This equation confirms that $\phi_i$ is an eigenvector of $\hat A$ that corresponds to eigenvalue $\lambda_i$ of both $\tilde A$ and $\hat A$. -``` +where $w_{t+1}$ and $v_t$ are standard normal random vectors. -Also see {cite}`DDSE_book` (p. 238) +The three columns of $G$ are shaped like a **level** factor that moves all quantiles together, a **slope** factor that moves low and high quantiles in opposite directions, and a **curvature** factor that moves the tails relative to the middle. +Here is code that simulates the data. -### Decoder of $\check b$ as a linear projection +```{code-cell} ipython3 +def simulate_data(m=300, T=134, D_diag=(0.95, 0.80, 0.50), σ_v=0.5, seed=1234): + """ + Simulate an m x T data matrix driven by three latent AR(1) factors. + """ + rng = np.random.default_rng(seed) + q = np.linspace(0, 1, m) + G = np.column_stack([np.ones(m), # level + 2 * (q - 0.5), # slope + 6 * (q - 0.5)**2 - 0.5]) # curvature + d = np.array(D_diag) + k = len(d) + z = np.zeros((k, T)) + z[:, 0] = rng.normal(size=k) / np.sqrt(1 - d**2) + for t in range(T - 1): + z[:, t+1] = d * z[:, t] + rng.normal(size=k) + X_tilde = G @ z + σ_v * rng.normal(size=(m, T)) + return X_tilde, G, z, q +X_tilde, G, z, q = simulate_data() +X_tilde = X_tilde - X_tilde.mean(axis=1, keepdims=True) # demean each variable +X_tilde.shape +``` +The next function implements the DMD algorithm. + +```{code-cell} ipython3 +def dmd(X_tilde, r): + """ + Rank-r DMD of the m x (n+1) data matrix X_tilde. + """ + X, X_prime = X_tilde[:, :-1], X_tilde[:, 1:] + U, σ, Vt = np.linalg.svd(X, full_matrices=False) + U_r, σ_r, V_r = U[:, :r], σ[:r], Vt[:r, :].T + + B = X_prime @ V_r / σ_r # X' V_r Σ_r^{-1} + A_tilde = U_r.T @ B # r x r VAR for principal components + λ, W = np.linalg.eig(A_tilde) + order = np.argsort(-np.abs(λ)) # most persistent modes first + λ, W = λ[order], W[:, order] + + Φ = B @ W # DMD modes (right eigenvectors) + Ψ = np.linalg.solve(W * λ, U_r.T) # (W Λ)^{-1} U_r^T (left eigenvectors) + A_hat = B @ U_r.T # rank-r estimator of A + a_hat = X_prime - A_hat @ X # residuals + Ω_hat = a_hat @ a_hat.T / X.shape[1] + + return dict(λ=λ, Φ=Φ, Ψ=Ψ, W=W, A_tilde=A_tilde, A_hat=A_hat, + Ω_hat=Ω_hat, σ=σ, U_r=U_r) +``` +To choose $r$, we plot the largest singular values of $X$. +```{code-cell} ipython3 +σ = np.linalg.svd(X_tilde[:, :-1], compute_uv=False) +fig, ax = plt.subplots() +ax.plot(np.arange(1, 21), σ[:20], 'o') +ax.set_xlabel('index $j$') +ax.set_ylabel(r'singular value $\sigma_j$') +plt.show() +``` -From eigendecomposition {eq}`eq:APhiLambda` we can represent $\hat A$ as +The singular values fall rapidly for $j = 1, 2, 3$ and then level off, so the elbow rule tells us to set $r = 3$. -$$ -\hat A = \Phi \Lambda \Phi^+ . -$$ (eq:Aform12) +```{code-cell} ipython3 +r = 3 +res = dmd(X_tilde, r) +λ, Φ, Ψ, A_hat = res['λ'], res['Φ'], res['Ψ'], res['A_hat'] +print("DMD eigenvalues: ", np.round(λ, 3)) +print("true eigenvalues: ", np.array([0.95, 0.80, 0.50])) +``` -From formula {eq}`eq:Aform12` we can deduce dynamics of the $p \times 1$ vector $\check b_t$: +The two most persistent eigenvalues are estimated well. -$$ -\check b_{t+1} = \Lambda \check b_t -$$ +The least persistent one is biased toward zero. -where +That is because the data are noisy measurements of the factors, so $\{X_t\}$ is not exactly a first-order VAR; errors in the principal components $f_t$ attenuate the estimated coefficients, most visibly for the least persistent mode. -$$ -\check b_t = \Phi^+ X_t -$$ (eq:decoder102) +Next we verify {prf:ref}`prop-dmd-eigenvectors` and {prf:ref}`prop-dmd-left`, and check whether $\Phi \Lambda \Phi^+$ equals $\hat A_r$. +```{code-cell} ipython3 +Φ_pinv = np.linalg.pinv(Φ) +rel = lambda M: np.linalg.norm(M) / np.linalg.norm(A_hat) -Since the $m \times p$ matrix $\Phi$ has $p$ linearly independent columns, the generalized inverse of $\Phi$ is +print("max |A_hat Φ - Φ Λ| =", np.abs(A_hat @ Φ - Φ * λ).max()) +print("max |Ψ Φ - I| =", np.abs(Ψ @ Φ - np.eye(r)).max()) +print("max |Ψ A_hat - Λ Ψ| =", np.abs(Ψ @ A_hat - λ[:, None] * Ψ).max()) +print("relative error, Φ Λ Ψ =", rel(A_hat - (Φ * λ) @ Ψ)) +print("relative error, Φ Λ Φ^+ =", rel(A_hat - (Φ * λ) @ Φ_pinv)) +``` -$$ -\Phi^{+} = (\Phi^\top \Phi)^{-1} \Phi^\top -$$ +The first four numbers are zero up to rounding error, while $\Phi \Lambda \Phi^+$ misses $\hat A_r$ by several percent. -and so +Now let's look at the modes. -$$ -\check b = (\Phi^\top \Phi)^{-1} \Phi^\top X -$$ (eq:checkbform) +Modes are determined only up to scale, so we normalize each mode to have maximum absolute value one and choose its sign to line up with the corresponding column of $G$, which we normalize in the same way. -The $p \times n$ matrix $\check b$ is recognizable as a matrix of least squares regression coefficients of the $m \times n$ matrix -$X$ on the $m \times p$ matrix $\Phi$ and consequently +```{code-cell} ipython3 +def normalize(v): + return v / np.abs(v).max() -$$ -\check X = \Phi \check b -$$ (eq:Xcheck_) +labels = ['level', 'slope', 'curvature'] +Φ_real = np.real_if_close(Φ) -is an $m \times n$ matrix of least squares projections of $X$ on $\Phi$. +fig, axes = plt.subplots(1, 3, figsize=(12, 3.5)) +for i, ax in enumerate(axes): + g_i = normalize(G[:, i]) + φ_i = normalize(Φ_real[:, i]) + φ_i *= np.sign(φ_i @ g_i) + ax.plot(q, φ_i, lw=2, label=rf'mode $\phi_{i+1}$') + ax.plot(q, g_i, 'k--', label=f'true {labels[i]} loading') + ax.set_ylim(-1.1, 1.1) + ax.set_xlabel('quantile $q$') + ax.legend() +plt.tight_layout() +plt.show() +``` - **Variance Decomposition of $X$** +Each mode is close to one of the true loading patterns, but not identical to it. -By virtue of the least-squares projection theory discussed in this quantecon lecture , we can represent $X$ as the sum of the projection $\check X$ of $X$ on $\Phi$ plus a matrix of errors. +Because of sampling error, the estimated eigenvectors mix the true factors a little, most visibly for the second mode. +Next we compare the modal coordinates $b_t = \Psi X_t$ with the true factors $z_t$, after standardizing both. -To verify this, note that the least squares projection $\check X$ is related to $X$ by +```{code-cell} ipython3 +b = np.real_if_close(Ψ @ X_tilde) +standardize = lambda x: (x - x.mean()) / x.std() +fig, axes = plt.subplots(3, 1, figsize=(10, 7), sharex=True) +for i, ax in enumerate(axes): + corr = np.corrcoef(b[i], z[i])[0, 1] + ax.plot(np.sign(corr) * standardize(b[i]), lw=2, label=rf'$b_{{{i+1},t}}$') + ax.plot(standardize(z[i]), 'k--', label=rf'$z_{{{i+1},t}}$') + ax.set_ylabel(f'corr = {abs(corr):.2f}') + ax.legend(loc='upper right') +axes[-1].set_xlabel('$t$') +plt.tight_layout() +plt.show() +``` -$$ -X = \check X + \epsilon -$$ +The DMD recovers the latent factors well. -or +Finally, we compare $j$-step-ahead forecasts computed in three ways: -$$ -X = \Phi \check b + \epsilon -$$ (eq:Xbcheck) +* directly from $\hat A_r^j X_t$ +* from the modal representation {eq}`eq:modalforecast` with $b_t = \Psi X_t$ +* from the tempting alternative $\Phi \Lambda^j \Phi^+ X_t$ -where $\epsilon$ is an $m \times n$ matrix of least squares errors satisfying the least squares orthogonality conditions $\epsilon^\top \Phi =0 $ or +```{code-cell} ipython3 +X = X_tilde[:, :-1] +for j in (1, 4, 12): + direct = np.linalg.matrix_power(A_hat, j) @ X + via_Ψ = np.real_if_close((Φ * λ**j) @ Ψ @ X) + via_pinv = np.real_if_close((Φ * λ**j) @ Φ_pinv @ X) + err = lambda F: np.linalg.norm(F - direct) / np.linalg.norm(direct) + print(f"j = {j:2d}: error with Ψ = {err(via_Ψ):.1e}, " + f"error with Φ^+ = {err(via_pinv):.3f}") +``` -$$ -(X - \Phi \check b)^\top \Phi = 0_{m \times p} -$$ (eq:orthls) +Forecasts built from $\Psi$ reproduce the VAR forecasts exactly, while those built from $\Phi^+$ do not. -Rearranging the orthogonality conditions {eq}`eq:orthls` gives $X^\top \Phi = \check b \Phi^\top \Phi$, which implies formula {eq}`eq:checkbform`. +In this example the forecast discrepancy from using $\Phi^+$ is small, but, as {ref}`var_dmd_ex1` shows, it grows with the amount of noise in the data. +## Source for some Python code +You can find a Python implementation of DMD in the [PyDMD](https://github.com/PyDMD/PyDMD) package. +## Exercises +```{exercise} +:label: var_dmd_ex1 -### An Approximation +In this lecture we saw that $\Phi^+ = \Psi$ if and only if the column spaces of $\Phi$ and $U_r$ are the same. +1. Simulate data with `simulate_data(σ_v=0.0)`, so that the data are exact linear combinations of three factors, and compute a rank-3 DMD. + Verify numerically that $\Phi^+ = \Psi$ and that $\Phi \Lambda \Phi^+ = \hat A_r$. +1. Explain why the column spaces of $\Phi$ and $U_r$ coincide in this case. +1. Repeat part 1 with $\sigma_v = 0.1, 0.5, 1.0$ and report how far $\Phi \Lambda \Phi^+$ is from $\hat A_r$. +``` +```{solution-start} var_dmd_ex1 +:class: dropdown +``` -We now describe a way to approximate the $p \times 1$ vector $\check b_t$ instead of using formula {eq}`eq:decoder102`. +Here is one solution. + +```{code-cell} ipython3 +for σ_v in (0.0, 0.1, 0.5, 1.0): + X_sim = simulate_data(σ_v=σ_v)[0] + X_sim = X_sim - X_sim.mean(axis=1, keepdims=True) + out = dmd(X_sim, 3) + Φ_s, Ψ_s, λ_s, A_s = out['Φ'], out['Ψ'], out['λ'], out['A_hat'] + Φ_s_pinv = np.linalg.pinv(Φ_s) + gap_left = np.linalg.norm(Φ_s_pinv - Ψ_s) / np.linalg.norm(Ψ_s) + gap_A = np.linalg.norm(A_s - (Φ_s * λ_s) @ Φ_s_pinv) / np.linalg.norm(A_s) + print(f"σ_v = {σ_v:.1f}: |Φ^+ - Ψ|/|Ψ| = {gap_left:.1e}, " + f"|A_hat - Φ Λ Φ^+|/|A_hat| = {gap_A:.1e}") +``` -In particular, the following argument adapted from {cite}`DDSE_book` (page 240) provides a computationally efficient way to approximate $\check b_t$. +When $\sigma_v = 0$, $X_t = G z_t$ for every $t$, so every column of both $X$ and $X'$ lies in the three-dimensional column space of $G$. -For convenience, we'll apply the method at time $t=1$. +Because $X$ has rank three, the column space of $U_3$ equals the column space of $G$. +The columns of $\Phi = X' V_3 \Sigma_3^{-1} W$ are linear combinations of columns of $X'$, so they also lie in the column space of $G$. +Since $\Phi$ has three linearly independent columns, its column space is the column space of $G$, which is the column space of $U_3$. -For $t=1$, from equation {eq}`eq:Xbcheck` we have +Then the rows of $\Psi$ lie in the space spanned by the rows of $\Phi^*$. -$$ - \check X_1 = \Phi \check b_1 -$$ (eq:X1proj) +So the rows of $\Psi - \Phi^+$ can be written as $C \Phi^*$ for some matrix $C$. -where $\check b_1$ is a $p \times 1$ vector. +Because $(\Psi - \Phi^+) \Phi = I - I = 0$, we get $C \Phi^* \Phi = 0$, and since $\Phi^* \Phi$ is invertible, $C = 0$. -Recall from representation 1 above that $X_1 = U \tilde b_1$, where $\tilde b_1$ is a time $1$ basis vector for representation 1 and $U$ is from the full SVD $X = U \Sigma V^\top$. +With measurement noise, $X'$ has components outside the column space of $U_3$, the two left inverses differ, and the gap between $\Phi \Lambda \Phi^+$ and $\hat A_r$ grows with $\sigma_v$. -It then follows from equation {eq}`eq:Xbcheck` that +```{solution-end} +``` - -$$ - U \tilde b_1 = X' \tilde V \tilde \Sigma^{-1} \tilde W \check b_1 + \epsilon_1 -$$ +```{exercise} +:label: var_dmd_ex2 -where $\epsilon_1$ is a least-squares error vector from equation {eq}`eq:Xbcheck`. +This exercise compares the DMD estimator $\hat A_r$ with the reduced-rank regression estimator of {cite}`anderson1951`. -It follows that +The reduced-rank regression estimator is the rank-$r$ matrix that minimizes $\| X' - \check A X \|_F$. -$$ -\tilde b_1 = U^\top X' V \tilde \Sigma^{-1} \tilde W \check b_1 + U^\top \epsilon_1 -$$ +It can be computed as $\hat A^{RRR}_r = P_r \hat A$, where $\hat A = X' X^+$ is the minimum-norm least squares estimator and $P_r = Q_r Q_r^\top$, where the columns of $Q_r$ are the first $r$ left singular vectors of the fitted values $\hat A X$. +Using data from `simulate_data()`, estimate both rank-3 estimators on the first 100 periods, then -Replacing the error term $U^\top \epsilon_1$ by zero, and replacing $U$ from a **full** SVD of $X$ with $\tilde U$ from a **reduced** SVD, we obtain an approximation $\hat b_1$ to $\tilde b_1$: +1. compute the relative in-sample residual norm $\| X' - \check A X \|_F / \| X' \|_F$ for each estimator +1. compute the same statistic out of sample, using one-step-ahead forecasts for the remaining periods +1. repeat for several random seeds and comment on the results +``` +```{solution-start} var_dmd_ex2 +:class: dropdown +``` +Here is one solution. + +```{code-cell} ipython3 +def rrr(X, X_prime, r): + "Reduced-rank regression estimator of rank r." + A_full = X_prime @ np.linalg.pinv(X) + Q = np.linalg.svd(A_full @ X, full_matrices=False)[0][:, :r] + return Q @ Q.T @ A_full + +def rel_resid(A, X, X_prime): + return np.linalg.norm(X_prime - A @ X) / np.linalg.norm(X_prime) + +n_est = 100 +print("seed DMD in DMD out RRR in RRR out") +for seed in range(5): + X_sim = simulate_data(seed=seed)[0] + X_sim = X_sim - X_sim.mean(axis=1, keepdims=True) + X_in, Xp_in = X_sim[:, :n_est-1], X_sim[:, 1:n_est] + X_out, Xp_out = X_sim[:, n_est-1:-1], X_sim[:, n_est:] + + A_dmd = dmd(X_sim[:, :n_est], 3)['A_hat'] + A_rrr = rrr(X_in, Xp_in, 3) + print(f"{seed:4d} {rel_resid(A_dmd, X_in, Xp_in):.3f} " + f"{rel_resid(A_dmd, X_out, Xp_out):.3f} " + f"{rel_resid(A_rrr, X_in, Xp_in):.3f} " + f"{rel_resid(A_rrr, X_out, Xp_out):.3f}") +``` -$$ - \hat b_1 = \tilde U^\top X' \tilde V \tilde \Sigma^{-1} \tilde W \check b_1 -$$ +By construction, reduced-rank regression has the smaller in-sample residuals. -Recall that from equation {eq}`eq:tildeAverify`, $ \tilde A = \tilde U^\top X' \tilde V \tilde \Sigma^{-1}$. +But in these samples it forecasts worse out of sample than the DMD estimator. -It then follows that - -$$ - \hat b_1 = \tilde A \tilde W \check b_1 -$$ +Reduced-rank regression starts from $\hat A = X' V \Sigma^{-1} U^\top$, which uses **all** singular values of $X$, including tiny ones that amplify noise, and then chooses the three directions that best fit $X'$ in sample. -and therefore, by the eigendecomposition {eq}`eq:tildeAeigenred` of $\tilde A$, we have +The DMD estimator instead discards the small singular values before it regresses, which acts as a form of regularization. -$$ - \hat b_1 = \tilde W \Lambda \check b_1 -$$ +```{solution-end} +``` -Consequently, - -$$ - \hat b_1 = ( \tilde W \Lambda)^{-1} \tilde b_1 -$$ +```{exercise} +:label: var_dmd_ex3 -or +Using the rank-3 DMD of the simulated data from the lecture, +1. prove formula {eq}`eq:Ahatpower` and use it to derive formula {eq}`eq:longrun` +1. verify formula {eq}`eq:longrun` numerically by comparing it with a direct computation of $(I - \hat A_r)^{-1}$ +1. treating the first component of $X_t$ as the growth-rate variable, compute $g = e_1 \hat A_r (I - \hat A_r)^{-1}$ both directly and with the modal formula $g = e_1 \Phi \Lambda (I - \Lambda)^{-1} \Psi$, and report the contribution of each mode to $g$ +1. compute the modal formula again with $\Phi^+$ in place of $\Psi$ and report how much the answer changes +``` -$$ - \hat b_1 = ( \tilde W \Lambda)^{-1} \tilde U^\top X_1 , -$$ (eq:beqnsmall) +```{solution-start} var_dmd_ex3 +:class: dropdown +``` +For part 1, {prf:ref}`prop-dmd-left` gives $\hat A_r = \Phi \Lambda \Psi$ and $\Psi \Phi = I$. +So $\hat A_r^2 = \Phi \Lambda (\Psi \Phi) \Lambda \Psi = \Phi \Lambda^2 \Psi$, and induction gives $\hat A_r^j = \Phi \Lambda^j \Psi$ for $j \geq 1$. -which is a computationally efficient approximation to the following instance of equation {eq}`eq:decoder102` for the initial vector $\check b_1$: +When $|\lambda_i| < 1$ for all $i$, $$ - \check b_1= \Phi^{+} X_1 -$$ (eq:bphieqn) - - -(To highlight that {eq}`eq:beqnsmall` is an approximation, users of DMD sometimes call components of basis vector $\check b_t = \Phi^+ X_t $ the **exact** DMD modes and components of $\hat b_t = ( \tilde W \Lambda)^{-1} \tilde U^\top X_t$ the **approximate** modes.) - -Conditional on $X_t$, we can compute a decoded $\check X_{t+j}, j = 1, 2, \ldots $ from the exact modes via - +\sum_{j=1}^\infty \hat A_r^j = \Phi \left( \sum_{j=1}^\infty \Lambda^j \right) \Psi = \Phi \Lambda (I - \Lambda)^{-1} \Psi $$ -\check X_{t+j} = \Phi \Lambda^j \Phi^{+} X_t -$$ (eq:checkXevoln) - - -or use compute a decoded $\hat X_{t+j}$ from approximate modes via -$$ - \hat X_{t+j} = \Phi \Lambda^j (\tilde W \Lambda)^{-1} \tilde U^\top X_t . -$$ (eq:checkXevoln2) +and since all eigenvalues of $\hat A_r$ are inside the unit circle, $(I - \hat A_r)^{-1} = I + \sum_{j=1}^\infty \hat A_r^j$. -We can then use a decoded $\check X_{t+j}$ or $\hat X_{t+j}$ to forecast $X_{t+j}$. +Here is code for parts 2--4. +```{code-cell} ipython3 +m = A_hat.shape[0] +I = np.eye(m) +weights = λ / (1 - λ) +LR_direct = np.linalg.inv(I - A_hat) +LR_modal = np.real_if_close(I + (Φ * weights) @ Ψ) +print("max |direct - modal| for (I - A_hat)^{-1}:", + np.abs(LR_direct - LR_modal).max()) -### Using Fewer Modes +e1 = np.zeros(m) +e1[0] = 1.0 +g_direct = e1 @ A_hat @ LR_direct +g_modal = np.real_if_close(e1 @ (Φ * weights) @ Ψ) +print("max |g_direct - g_modal|:", np.abs(g_direct - g_modal).max()) -In applications, we'll actually use only a few modes, often three or less. +print("\nlong-run weights λ/(1-λ):", np.round(weights.real, 2)) +for i in range(r): + g_i = np.real_if_close(Φ[0, i] * weights[i] * Ψ[i, :]) + share = np.linalg.norm(g_i) / np.linalg.norm(g_direct) + print(f"mode {i+1}: |contribution to g| / |g| = {share:.2f}") -Some of the preceding formulas assume that we have retained all $p$ modes associated with singular values of $X$. - -We can adjust our formulas to describe a situation in which we instead retain only -the $r < p$ largest singular values. - -In that case, we simply replace $\tilde \Sigma$ with the appropriate $r\times r$ matrix of singular values, $\tilde U$ with the $m \times r$ matrix whose columns correspond to the $r$ largest singular values, and $\tilde V$ with the $n \times r$ matrix whose columns correspond to the $r$ largest singular values. - -Counterparts of all of the salient formulas above then apply. +g_pinv = np.real_if_close(e1 @ (Φ * weights) @ Φ_pinv) +print("\nrelative error in g using Φ^+:", + np.linalg.norm(g_pinv - g_direct) / np.linalg.norm(g_direct)) +``` +The modal formula reproduces the direct computation up to rounding error, and requires only $r \times r$ operations beyond forming $\Phi$ and $\Psi$. +The most persistent mode, with the largest weight $\lambda_i/(1-\lambda_i)$, contributes the most to $g$. -## Source for Some Python Code +The ratios need not sum to one, and one can exceed one, because contributions of different modes partially offset each other. -You can find a Python implementation of DMD here: +Replacing $\Psi$ by $\Phi^+$ gives a different, incorrect answer. -https://mathlab.sissa.it/pydmd +```{solution-end} +``` From 051b9ba9948d70d1d6eda472874792b570046a47 Mon Sep 17 00:00:00 2001 From: John Stachurski Date: Tue, 15 Sep 2026 18:09:44 +1000 Subject: [PATCH 2/2] [var_dmd, svd_intro, measurement_models] Review fixes: eigenvalue printouts, docstring timing, heading case MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - measurement_models: print V2 eigenvalues in scientific notation so the tiny second eigenvalue and the third (= ε) are visible, as the text claims - measurement_models: relabel the Friedman step in the meas_ex1 solution as the second claim; align steady_state_kalman docstring with the w_t timing convention - svd_intro: sentence-case section headings per the style manual - var_dmd: i.i.d. -> IID - quant-econ.bib: add trailing newline Co-Authored-By: Claude Opus 5 (1M context) --- lectures/_static/quant-econ.bib | 2 +- lectures/measurement_models.md | 9 +++++---- lectures/svd_intro.md | 20 ++++++++++---------- lectures/var_dmd.md | 2 +- 4 files changed, 17 insertions(+), 16 deletions(-) diff --git a/lectures/_static/quant-econ.bib b/lectures/_static/quant-econ.bib index 86bf3c246..824a1d029 100644 --- a/lectures/_static/quant-econ.bib +++ b/lectures/_static/quant-econ.bib @@ -5404,4 +5404,4 @@ @article{Geweke1982 number = {378}, pages = {304--313}, year = {1982} -} \ No newline at end of file +} diff --git a/lectures/measurement_models.md b/lectures/measurement_models.md index 1746a53b5..ab8bb0f5a 100644 --- a/lectures/measurement_models.md +++ b/lectures/measurement_models.md @@ -503,7 +503,7 @@ innovation covariance $V$ def steady_state_kalman(A, C_obs, Q, R, W=None, tol=1e-13, max_iter=200_000): """ Solve steady-state Kalman equations for - x_{t+1} = A x_t + w_{t+1} + x_{t+1} = A x_t + w_t y_t = C_obs x_t + v_t with cov(w)=Q, cov(v)=R, cov(w,v)=W. """ @@ -1233,7 +1233,8 @@ display(Latex(df_to_latex_matrix(df_v2))) ```{code-cell} ipython3 print('Eigenvalues of covariance matrix:') -print(np.sort(np.linalg.eigvalsh(V2))[::-1].round(4)) +print(np.array2string(np.sort(np.linalg.eigvalsh(V2))[::-1], + formatter={'float': '{:.2e}'.format})) ``` One eigenvalue is large. @@ -1576,7 +1577,7 @@ because $(f-1)\beta = 1 - \beta$. Similarly, $y_{nt} - y_{n,t-1} = (f-1)\beta \theta_{t-1} + \theta_t - \theta_{t-1} = \theta_t - \beta\theta_{t-1}$. -For the third claim, substitute $k_t = k_{t-1} + \beta\theta_{t-1}$ into {eq}`simple_crule` and {eq}`net_income`: +For the second claim, substitute $k_t = k_{t-1} + \beta\theta_{t-1}$ into {eq}`simple_crule` and {eq}`net_income`: $$ \begin{aligned} @@ -1704,7 +1705,7 @@ print("max |identity residual| in reported data:", for ε_R in (1e-2, 1e-4, 1e-6): _, _, V2_ε = steady_state_kalman(A, C, Q2, ε_R * np.eye(3)) print(f"ε = {ε_R:.0e}: eigenvalues of V2 = " - f"{np.sort(np.linalg.eigvalsh(V2_ε))[::-1]}") + f"{np.array2string(np.sort(np.linalg.eigvalsh(V2_ε))[::-1], formatter={'float': '{:.2e}'.format})}") ``` One eigenvalue of $V_2$ is large: it reflects news about the permanent shock $\theta_t$. diff --git a/lectures/svd_intro.md b/lectures/svd_intro.md index 4aca757af..fcb8ce829 100644 --- a/lectures/svd_intro.md +++ b/lectures/svd_intro.md @@ -29,7 +29,7 @@ In a sequel to this lecture, {doc}`VARs and DMDs `, we'll describe how Like principal components analysis, DMD can be thought of as a data-reduction procedure that represents salient patterns by projecting data onto a limited set of factors. -## The Setting +## The setting Let $X$ be an $m \times n$ matrix of rank $p$. @@ -59,7 +59,7 @@ In the $m \gg n$ case in which there are many more attributes $m$ than individua In the sequel {doc}`VARs and DMDs `, we'll again use a singular value decomposition, but now to construct a **dynamic mode decomposition** (DMD). -## Singular Value Decomposition +## Singular value decomposition A **singular value decomposition** of an $m \times n$ matrix $X$ of rank $p \leq \min(m,n)$ is @@ -115,7 +115,7 @@ Thus, We'll apply this circle of ideas in the sequel {doc}`VARs and DMDs ` when we study dynamic mode decompositions. -**Road Ahead** +**Road ahead** What we have described above is called a **full** SVD. @@ -125,7 +125,7 @@ Later we'll also describe an **economy** or **reduced** SVD. Before we study a **reduced** SVD we'll say a little more about properties of a **full** SVD. -## Four Fundamental Subspaces +## Four fundamental subspaces Let ${\mathcal C}$ denote a column space, ${\mathcal N}$ denote a null space, and ${\mathcal R}$ denote a row space. @@ -331,7 +331,7 @@ print("col_space ⟂ left_null_space:", np.allclose(col_space.T @ left_null_spac print("row_space ⟂ null_space: ", np.allclose(row_space.T @ null_space, 0)) ``` -## Eckart-Young Theorem +## Eckart-Young theorem Suppose that we want to construct the best rank $r$ approximation of an $m \times n$ matrix $X$. @@ -387,7 +387,7 @@ of the approximation $\hat X_r$ is a key ingredient of a dynamic mode decomposit {ref}`svd_ex2` and {ref}`svd_ex3` explore these ideas. -## Full and Reduced SVD's +## Full and reduced SVD's Up to now we have described properties of a **full** SVD in which shapes of $U$, $\Sigma$, and $V$ are $\left(m, m\right)$, $\left(m, n\right)$, $\left(n, n\right)$, respectively. @@ -536,7 +536,7 @@ SShat=np.diag(Shat) np.allclose(X, Uhat@SShat@Vhat) ``` -## Polar Decomposition +## Polar decomposition A **reduced** singular value decomposition (SVD) of $X$ is related to a **polar decomposition** of $X$ @@ -566,7 +566,7 @@ Here * $S$ is an $m \times m$ **symmetric positive semidefinite** matrix * $Q$ is an $m \times n$ matrix with $Q Q^\top = U U^\top$; so when $p = m$, as for a short-fat $X$ with full row rank, $Q Q^\top = I$ and $Q$ has orthonormal rows -## Application: Principal Components Analysis (PCA) +## Application: principal components analysis (PCA) Let's begin with a case in which $n \gg m$, so that we have many more individuals $n$ than attributes $m$. @@ -659,7 +659,7 @@ $$ By the Eckart-Young theorem, one minus this ratio is also the squared relative Frobenius-norm error $||B - \hat B_r||_F^2 / ||B||_F^2$ of the best rank $r$ approximation $\hat B_r$ of $B$. -**Step 5: Create the Score Matrix:** +**Step 5: Create the score matrix:** The matrix of **principal components** (or **scores**) is @@ -717,7 +717,7 @@ is a vector of **loadings** of variables $X_i$ on the $k$th principal component, * $\sigma_k $ for each $k=1, \ldots, p$ is the strength of the $k$th **principal component**, where strength means contribution to the overall covariance of $X$: the $k$th principal component contributes $\sigma_k^2 / n$ to the sum of the sample variances of the $m$ variables -## PCA with Eigenvalues and Eigenvectors +## PCA with eigenvalues and eigenvectors We now use an eigen decomposition of a sample covariance matrix to do PCA. diff --git a/lectures/var_dmd.md b/lectures/var_dmd.md index 67196c145..f2f64d734 100644 --- a/lectures/var_dmd.md +++ b/lectures/var_dmd.md @@ -56,7 +56,7 @@ $$ X_{t+1} = A X_t + C \epsilon_{t+1}, \quad \epsilon_{t+1} \perp X_t $$ (eq:VARfirstorder) -where $\epsilon_{t+1}$ is the time $t+1$ component of a sequence of i.i.d. $m \times 1$ random vectors with mean vector zero and identity covariance matrix and where the $m \times 1$ vector $X_t$ is +where $\epsilon_{t+1}$ is the time $t+1$ component of a sequence of IID $m \times 1$ random vectors with mean vector zero and identity covariance matrix and where the $m \times 1$ vector $X_t$ is $$ X_t = \begin{bmatrix} X_{1,t} & X_{2,t} & \cdots & X_{m,t} \end{bmatrix}^\top