Merge pull request #8 from VFCI/7-max-share-derivations

Updated draft of derivations

vignettes/articles/derivations.Rmd

@@ -1,96 +1,216 @@
 ---
-title: "Derivations for Frequnecy Domain FEV"
+title: "Derivations for Max Share Identification Method"
+author: "Matthew DeHaven"
+date: "January 16, 2023"
 ---
 
-```{r, include = FALSE}
-knitr::opts_chunk$set(
-  collapse = TRUE,
-  comment = "#>"
-)
-```
+Article explaining the mathematics behind the two key `fevdid` functions:
 
+- `id_fevdfd()`
+- `id_fevdtd()`
 
-# Vector Autoregression (VAR)
+## Structural Vector Autoregression (SVAR)
 
-Assume a VAR(p) process with the following form:
+A SVAR(p) model with $p$ lags, for a vector of variables $y_t$,
+assumes the data generating process is accurately represented as,
+^[Notation and setup follows Chapter 1: Introduction from Kilian and Lutkepohl (2017)]
 
 $$
-A(L) X_t = \upsilon_t
+B_0 y_t = B_1 y_{t-1} + ... + B_p y_{t-p} + w_t
 $$
 
-where 
+where $w_t$ are mean zero structural shocks, with a serially uncorrelated diagonal covariance matrix of $\Sigma_w = E[w_t w_t]$,
 
-- $X_t$ is an $N \times 1$ vector of the data
+and $B_i$ are matrices of coefficients.
 
-- $A(L)$ is a matrix polynomial with the lag operator, $L$
+Empirically, however, only the following $A_i$ matrices of coefficients and $v_t$ vector of reduced form residuals are observed,
 
-- $\upsilon_t$ is a vector of residuals, i.e. reduced-form shocks
+$$
+y_t = \underbrace{B_0^{-1}B_1}_{A_1} y_{t-1} + ... + \underbrace{B_0^{-1}B_p}_{A_p} y_{t-p} + \underbrace{B_0^{-1} w_t}_{u_t} 
+\\
+y_t = A_1 y_{t-1} + ... + A_p y_{t-p} + u_t 
+$$
+
+#### Identification Problem
+
+We observe reduced form residuals, $v_t$, but want the structural shocks, $w_t$.
+
+$$
+u_t = B_0^{-1} w_t
+$$
+
+At this point, WLOG, assume that $\Sigma_w = I$.
 
-- $p$ represents the number of lags of the data included in the VAR
+This implies that $\Sigma_u = B_0^{-1} B_0^{-1^\prime}$.
 
-- $E[\upsilon_t \upsilon_t'] = \Sigma$, so $\Sigma$ is the variance-covariance matrix of the residuals
+In order to identify the structural shocks, we simply need to know the structural impact matrix, $B_0$, or equivalently, the inverse $B_0^{-1}$.
 
-with the lag matrix having the following structure:
+### Impulse Response Functions (IRF) and Forecast Errors
 
+It will be useful for the max share method to show the impulse response and forecast error definitions
+
+#### IRFs
+IRFs give the responses of each variable in $y_t$ for any horizon $i$ to a one-time impulse at time $t$ to each structural shock, $w_t$,
+$$
+\Theta_i \equiv \frac{\partial y_{t+i}}{\partial w_t} \hspace{1cm} i = 0, 1, 2, ..., H
+$$
+
+The IRF for the response of variable $k$ to an impulse for structural shock $j$ at horizon $i$ is denoted,
 $$
-A(L) \equiv \sum_{T=0}^{p} A_T L^T
+\theta_{kj,i} \equiv \frac{\partial y_{k, t+i}}{\partial w_{jt}}
+$$
+
+Let the impulse responses to the reduced form residuals be similarily denoted as,
+$$
+\Phi_i \equiv \frac{\partial y_{t+i}}{\partial u_t} \hspace{1cm} i = 0, 1, 2, ..., H
 \\
-A(0) = A_0 = I
+\phi_{kj,i} \equiv \frac{\partial y_{k, t+i}}{\partial u_{jt}}
 $$
 
-# Structural shocks
+And we can map between the two using the $B_0$ impact matrix,
+$$
+\Theta_i = \Phi_i B_0^{-1}
+\\
+\theta_{kj,i} = [\Phi_i B_0^{-1}]_{kj}
+$$
 
-Assume there is a linear mapping, $S$, between $v_t$ and some independent structural shocks, $\epsilon_t$.
+Letting $\Phi_{k*,i}$ denote the $k$th row of $\Phi_i$ and $B_{0,*j}^{-1}$ denote the $j$th column of $B_0^{-1}$, then
 
 $$
-\upsilon_t = S \epsilon_t
+\theta_{kj,i} = \Phi_{k*,i} B_{0,*j}^{-1}
 $$
 
-with $E[\epsilon_t \epsilon_t'] = I$.
+#### Forecast Errors
+
+Let $F_{t+h}$ denote $h$-step ahead forecast errors,
+$$
+F_{t+h} \equiv y_{t+h} - y_{t+h|t}
+$$
 
-It is convenient to represent $S$ as a rotation of the Choleskey decomposition, $\tilde{S}$, of $\Sigma$.
+Which, for a VAR process, can be denoted in terms of the redcued form IRFs, $\Phi_i$,
 $$
-S = \tilde{S}Q
+F_{t+h} = \sum_{i=0}^{h-1} \Phi_i u_{t+h-i} = \sum_{i=0}^{h-1} \Phi_i B_0^{-1} w_{t+h-i}
 $$
-where $Q$ is an orthonormal matrix ($Q^{-1} = Q'$) so that it represents a rotation.
 
-There are infinitely many possible rotations $Q$ to satisfy the linear mapping.
-Taking $Q = I$ uses the Choleskey decomposition for identification and relies upon the ordering of the variables.
-Sign restriction identification involes drawing many possible matrices $Q$ and keeping only the set that have the required signs on impact.
-The Angeletos, Collard, and Dellas (2020) paper and this package instead identify the $Q$ matrix that leads to one structural shock explaining the largest portion of the forecast error variance for a targetted variable, either in the time or frequency domain.
 
-# Moving Average (Wold) Representation
+#### Forecast Error Variance (FEV)
+The FEV, also known as the predicted mean squared error, is then calculated as
+$$
+FEV_h \equiv E[F_{t+h}F_{t+h}^\prime] 
+\\
+= \sum_{i=0}^{h-1} \Phi_i \Sigma_u \Phi_i^\prime
+= \sum_{i=0}^{h-1} \Phi_i B_0^{-1} B_0^{-1^\prime} \Phi_i^\prime 
+= \sum_{i=0}^{h-1} [\Phi_i B_0^{-1}] [\Phi_i B_0^{-1}]^\prime
+$$
 
-The moving average representation of the VAR is given by
+The FEV for impulse $j$ contribution to variable $k$ is then denoted as,
+$$
+FEV_{kj,h} =  \sum_{i=0}^{h-1}  [\Phi_i B_0^{-1}]_{jk} [\Phi_i B_0^{-1}]_{jk}^\prime
+\\
+= \sum_{i=0}^{h-1} \Phi_{k*,i} B_{0,*j}^{-1} [\Phi_{k*,i} B_{0,*j}^{-1}]^\prime
+\\
+= \sum_{i=0}^{h-1} \Phi_{k*,i} B_{0,*j}^{-1} B_{0,*j}^{-1\prime} \Phi_{k*,i}^\prime
+$$
 
+Since $\Phi_{k*,i} B_{0,*j}^{-1}$ is a single number (not a matrix), we can rearrange to
 $$
-X_t = B(L) \upsilon_t
+= \sum_{i=0}^{h-1} B_{0,*j}^{-1\prime} \Phi_{k*,i}^\prime \Phi_{k*,i} B_{0,*j}^{-1}
+\\
+= B_{0,*j}^{-1\prime} \sum_{i=0}^{h-1}\left[\Phi_{k*,i}^\prime \Phi_{k*,i}\right] B_{0,*j}^{-1}
 $$
 
-where $B(L) = A(L)^{-1} = \sum_{T=0}^{\infty}B_TL^T$ and is an infinite matrix polynomial.
+## Max Share Identification Method - Time Domain
 
-With this, it is easy to rewrite our model as
+Pick a structural shock $j$ by choosing the weights for $B_{0,*j}^{-1}$ that maximize the forecast error variance for the horizon in $[h^-, h^+]$
 
+Solve the maximization problem:
 $$
-X_t = C(L) Q \epsilon_t = \Gamma(L) \epsilon_t
+\max_{B_{0,*j}^{-1}} \sum_{h^-}^{h^+} \text{FEV}_{kj,h}
 $$
 
-where $C(L) = = \sum_{T=0}^{\infty}C_TL^T$ with $C_T \equiv B_T \tilde{S}$ and 
-$\Gamma(L) = = \sum_{T=0}^{\infty}\Gamma_TL^T$ with $\Gamma_T \equiv B_T \tilde{S} Q$.
+This is equivalent to solving
+$$
+\max_{B_{0,*j}^{-1}} B_{0,*j}^{-1\prime} \sum_{h^-}^{h^+} \sum_{i=0}^{h-1}\left[\Phi_{k*,i}^\prime \Phi_{k*,i}\right] B_{0,*j}^{-1}
+$$
 
-# Targetting the Time Domain
+Which is maximized when $B_{0,*j}^{-1}$ is the eigenvector associated with the largest eigenvalue of the matrix $\sum_{h^-}^{h^+} \sum_{i=0}^{h-1}\left[\Phi_{k*,i}^\prime \Phi_{k*,i}\right]$.
 
-# Targetting the Frequency Domain
+**There's a piece missing here, which is that $B_{0,*j}^{-1}B_{0,*j}^{-1}$ can be decomposed into a Choleskey matrix and rotation matrix Q, then the maximization problem becomes picking a column of Q.**
 
-# Representing a VAR(p) as a VAR(1)
+## Max Share Identification Method - Frequency Domain
 
+Pick a structural shock $j$ by choosing the weights for $B_{0,*j}^{-1}$ that maximize the forecast error variance for the frequencies in $[\omega^-, \omega^+]$
+
+Solve the maximization problem:
+$$
+\max_{B_{0,*j}^{-1}} \int_{\omega^-}^{\omega^+} \text{FEV}_{kj,\omega} d\omega
+$$
+
+This is equivalent to solving
+$$
+\max_{B_{0,*j}^{-1}} B_{0,*j}^{-1\prime} \int_{\omega^-}^{\omega^+} \Phi_{k*,\omega}^\prime \Phi_{k*,\omega} \ d\omega \ B_{0,*j}^{-1}
+$$
+
+where $\Phi_{k*,\omega}$ is the impulse response function for a specific frequency $\omega$.  The next sections show how to calculate this value.
+
+### Solving for IRF in frequency space
+Note: One approximate solution is to calculate the time domain IRF out to long horizon (say 1000 periods), then take the fourier transformation, keep only the frequencies you want, square the contributions, and maximize the real value portion.
+
+#### Representing VAR(p) as a VAR(1)
+
+First it's easiest to work with the VAR(1) process instead of the VAR(p)
+^[see page 25 of Killian and Lutkepohl]
+
+. . .
+
+$$
+Y_t = \upsilon + A Y_{t-1} + U_t
+$$
+
+. . .
+
+#### MA($\infty$) Representation
+Next, represent as a moving average.
+^[see page 26 of Killian and Lutkepohl]
+
+. . .
+
+$$
+y_t = C(L)^{-1} \upsilon + C(L)^{-1} u_t
+$$
+
+where $C(L)^{-1} \equiv \sum_{i=0}^{\infty} \phi_i L^i$
+
+. . .
+
+#### Calculating the spectrum for an MA($\infty$)
+Now, we can calculate the spectrum using
+^[See Hamilton, page 152]
+
+. . .
+
+$$
+s_Y(\omega) = \frac{1}{2\pi} \sigma^2 C(e^{-i\omega})C(e^{-i\omega})
+$$
+
+The key here is that we end up getting to,
+$$
+C(e^{-i\omega}) = (I - A_1 e^{-i\omega})^{-1}
+$$
+
+which is essentially solving out for the entire path of IRFs by inverting, and making the Fourier transformation at the same time.
+
+. . .
 
 ---
 
-# Sources
+# Resources
+
+Hamilton, James D. *Time Series Analsysis*. Chapter 6: Spectral Analysis. Pg 152 - 172.  In particular, page 154 for deriving the population spectrum for a MA(infinity) process.
 
-- Hamilton, James D. *Time Series Analsysis*. Chapter 6: Spectral Analysis. Pg 152 - 172.  In particular, page 154 for deriving the population spectrum for a MA(infinity) process.
+Kilian, Lutz and Lutekepohl Helmut. *Structural Vector Autoregressive Analysis*. (2017). In particular,
 
-- Kilian, Lutz and Lutekepohl Helmut. *Structural Vector Autoregressive Analysis*. Chapter 2: Vector Autoregressive Models. Pg 19 - .  In particular, page 25 for representing a VAR(p) process as a VAR(1) process.
+- *Chapter 1: Introduction* for SVAR notation.
+- *Chapter 2: Vector Autoregressive Models*. Pg 19 - .  Specifically, page 25 for representing a VAR(p) process as a VAR(1) process.
 
-- Angeletos, George-Marios, Fabrice Collard, and Harris Dellas. (2020). "Businsess Cycle Anatomy". American Economic Review, 110(10): 3030 - 3070.  In particular, pages 3036 - 3037 for derivations.
+Angeletos, George-Marios, Fabrice Collard, and Harris Dellas. (2020). "Businsess Cycle Anatomy". American Economic Review, 110(10): 3030 - 3070.  In particular, pages 3036 - 3037 for derivations.