Econometric issues with Laubach and Williams' estimates of the natural rate of interest

Econometric issues with Laubach and Williams' estimates of the natural rate of interest footnote0

\thispagestyle{empty}

\ifthenelse{\equal{1}{1}}{ \secondtitle[2.166][\SecondTitlespace]{Econometric issues with Laubach and Williams' \\ estimates of the natural rate of interest}

}

\ifthenelse{\equal{0}{1}}

\let\oldref\ref \setcounter{page}{1} \setstretch{1.23}

\stepcounter{section} \Osubsection*{A.\arabic{section}. {Introduction }} \addcontentsline{toc}{subsection}{\currentname}

Since the global financial crisis, nominal interest rates have declined substantially to levels last witnessed in the early 1940s following the Great Depression. The academic as well as policy literature has attributed this decline in nominal interest rates to a decline in the natural rate of interest; namely, the rate of interest consistent with employment at full capacity and inflation at its target. In this literature, \citeauthor*{ holston.etal:2017}' (holston.etal:2017) estimates of the natural rate have become particularly influential and are widely regarded as a benchmark. The Federal Reserve Bank of New York (FRBNY) maintains an entire website dedicated to providing updates to \cites{holston.etal:2017} estimates of the natural rate, not only for the United States (U.S.), but also for the Euro Area, Canada and the United Kingdom (U.K.) (see \url{https://www.newyorkfed.org/research/policy/rstar}).

In \cites{holston.etal:2017} model, the natural rate of interest is defined as the sum of trend growth of output $g_{t}$ and `other factor' $ z_{t} $. This `other factor' $z_{t}$ is meant to capture various underlying structural factors such as savings/investment imbalances, demographic changes, and fiscal imbalances that influence the natural rate, but which are not captured by trend growth $g_{t}$. In (ref) below, I show filtered (as well as smoothed) estimates of \cites{holston.etal:2017} `other factor' $z_{t}$.\footnote{holston.etal:2017 do not show a plot of `other factor' $z_{t}$ on the FRBNY website (as of 22$^{nd}$ of June, 2020).}

\vsp[-2]

The dashed lines in (ref) show estimates obtained with data ending in 2017:Q1, while the solid lines are estimates based on data extended to 2019:Q2. The strong and persistent downward trending behaviour of `other factor' $z_{t}$ is striking from (ref), particularly from 2012:Q1 onwards. The two (black) dashed vertical lines mark the periods 2012:Q1 and 2015:Q4. In 2015:Q4, the Federal Reserve started the tightening cycle and raised nominal interest rates by 25 basis points. In 2012:Q1, real rates began to rise due to a (mild) deterioration in inflation expectations.\footnote{ See panel (a) of (ref), which shows plots of the federal funds rate, the real interest rate, as well as inflation and inflation expectations.} Both led to an increase in the real rate. Yet, \cites{holston.etal:2017} estimates of `other factor' $z_{t}$ declined by about 50 basis points from 2012:Q1 to 2015:Q4, and then another 50 basis points from 2015:Q4 to 2019:Q2, reaching a value of $-1.58$ in 2019:Q2. Because $z_{t}$ evolves as a driftless random walk in the model, the only parameter that `controls' the influence of $z_{t}$ on the natural rate is the `signal-to-noise ratio' $\lambda _{z}$.\footnote{ This description is somewhat imprecise to avoid cumbersome language. Since $ z_{t}$ evolves as $z_{t}=z_{t-1}+\sigma _{z}\epsilon _{t}$, with $\epsilon _{t}$ being standard normal, it is the standard deviation $\sigma _{z}$ that is the only parameter that influences the evolution of $z_{t}$. However, holston.etal:2017 determine $\sigma _{z}$ indirectly through the ` signal-to-noise ratio' $\lambda _{z}$, so it is the size of $\lambda _{z}$ that matters for the evolution of $z_{t}$.} Thus, how exactly this parameter is estimated is of fundamental importance for the determination of the natural rate of interest.

In this paper, I show that \cites{holston.etal:2017} implementation of \cites{stock.watson:1998} Median Unbiased Estimation (MUE) is unsound. It cannot recover the ratios of interest $\lambda _{g}=\sigma _{g}/\sigma _{y^{\ast }}$ and $\lambda _{z}=a_{r}\sigma _{z}/\sigma _{\tilde{y}}$ from Stages 1 and 2 of their three stage procedure needed for the estimation of the full structural model. The implementation of MUE of $\lambda _{z}$ in Stage 2 is particularly problematic, as \cites{holston.etal:2017} procedure is based on an `unnecessarily' misspecified Stage 2 model. This misspecified Stage 2 model not only fails to identify the ratio of interest $ \lambda _{z}=a_{r}\sigma _{z}/\sigma _{\tilde{y}}$, but moreover, due to the way holston.etal:2017 implement MUE in Stage 2, leads to spuriously large and excessively amplified estimates of $\lambda _{z}$. Since the magnitude of $\lambda _{z}$ determines and drives the downward trending behaviour of `other factor' $z_{t}$, this misspecification is consequential. Correcting their Stage 2 model and the MUE\ implementation results in a substantial quantitative reduction in the point estimate of $ \lambda _{z}$, and hence also $\sigma _{z}$. For instance, using data ending in 2017:Q1, \cites{holston.etal:2017} estimate of $\lambda _{z}$ is $ 0.030217 $ and yields an implied value of $0.150021$ for $\sigma _{z}$. After the correction, $\lambda _{z}$ is estimated to be $0.000754$ with an implied value for $\sigma _{z}$ of $0.003746$.\footnote{ These are my replicated estimates using data up to 2017:Q1, but they are effectively identical to those listed in Table 1, column 1 for the U.S. on page S60 in holston.etal:2017.} The resulting filtered (and smoothed) estimates of $z_{t}$ are markedly different, with the one from the correct Stage 2 implementation not only being very close to zero, but also highly insignificant statistically. The $p-$values corresponding to the structural break statistics from which $\lambda _{z}$ is estimated are of an order of magnitude of 0.5. These results highlight that there is no evidence of `other factor' $z_{t}$ in this model. The large and persistent downward trend in \cites{holston.etal:2017} estimates thus appears to be spurious.

In \Sref{sec:S2}, I outline in detail the Stage 2 model and the MUE procedure that holston.etal:2017 implement to estimate $\lambda _{z}$ . I show that their Stage 2 model is misspecified and that due to this, their MUE procedure cannot identify the ratio of interest $a_{r}\sigma _{z}/\sigma _{\tilde{y}}$ from $\lambda _{z}$. Instead, it recovers $\lambda _{z}=a_{r}\sigma _{z}/(\sigma _{\tilde{y}}+0.5a_{g}\sigma _{g})$ if $ (a_{g}+4a_{r})=0$. If $(a_{g}+4a_{r})\neq 0$, then additional parameters enter the denominator of $\lambda _{z}$, making it more intricate to recover $\sigma _{z}$ from $\lambda _{z}$, as it will be necessary to make additional assumptions about the time series properties of the nominal interest rate which is not explicitly modelled by holston.etal:2017, but rather added as an exogenous variable. The terms $a_{r}$ and $a_{g}$ are the parameters on the lagged real interest rate and lagged trend growth in the Stage 2 model of the output gap equation (see \Sref{sec:S2} for more details). In the full model, these are restricted so that $a_{g}=-4a_{r}$. In their specification of the Stage 2 model, holston.etal:2017 do not impose this restriction. Moreover, they include only one lag of trend growth $g_{t}$ in the output gap equation and, curiously, further add an intercept term to the specification that is not present in the full model (see equation (\oldref{S2:ytilde})). Since \cites{stock.watson:1998} MUE relies upon chow:1960 type structural break tests to estimate $\lambda _{z}$, these differences in the output gap specification lead to substantially larger $F$ statistics (see (ref) for a visual presentation) and therefore estimates of $\lambda _{z}$. To demonstrate that their misspecified Stage 2 model and MUE procedure leads to spurious and excessively large estimates of $\lambda _{z}$ when the true value is zero, I implement a simulation experiment in \Sref{sec:S2}. This simulation experiment shows that the mean estimate of $\lambda _{z}$ can be as high as $ 0.028842$, with a $45.7\%$ probability (relative frequency) of observing a value larger than estimated from the empirical data, when computed from simulated data which were generated from a model with the true $\lambda _{z}=0$. These simulation results are concerning, because they suggest that it is \cites{holston.etal:2017} MUE procedure itself that leads to the excessively large estimates of $\lambda _{z}$, rather than the size of the true $\lambda _{z}$ in the data.

Although \Sref{sec:S2} describes the core problem with \cites{holston.etal:2017} estimation procedure, there are other issues with the model and how it is estimated. Some of these are outlined in \Sref{sec:other}.\ For instance, \cites{holston.etal:2017} estimates of the natural rate, trend growth, `other factor' $z_{t}$ and the output gap are extremely sensitive to the starting date of the sample used to estimate the model. Estimating the model with data beginning in 1972:Q1 (or 1967:Q1) leads to negative estimates of the natural rate of interest toward the end of the sample period. These negative estimates are again driven purely by the exaggerated downward trending behaviour of `other factor' $z_{t}$ . The 1972:Q1 sample start was chosen to match the starting date used in the estimation of this model for the Euro Area. Out of the four countries that \cites{holston.etal:2017} model is fitted to, only the Euro Area estimates of the natural rate turn negative in 2013.\footnote{ Only the Euro Area estimates are based on a sample that starts in 1972:Q1, while the estimates for the U.K., Canada and the U.S. are based on samples starting in 1961:Q1.} The fact that it is also possible to generate such negative estimates of the natural rate from \cites{holston.etal:2017} model for the U.S. by simply adjusting the start of the estimation period to match that of the Euro Area data suggests that the model is far from robust, and therefore inappropriate for use in policy analysis. Furthermore, because Kalman Filtered estimates of the natural rate of interest will be moving averages of the observed variables that enter the state-space model, a circular or confounding relationship between the natural rate and the (nominal)\ policy rate will arise, because any central bank induced change in the policy target will be mechanically transferred to the natural rate via the Kalman Filtered estimate of the state vector. This makes it impossible to address `causal' questions regarding the relationship between natural rates and policy rates.

Median Unbiased Estimation is neither well known nor widely used at policy institutions. To give some background on the methodology, and to be able to understand why \cites{holston.etal:2017} implementation of MUE\ in Stage 2 is unsound, I provide a concise but important and informative overview of the methodology in \Sref{sec:MUE}.\ This section is essential for readers unfamiliar with the estimator. It reviews and summarises the conditions when it is likely to encounter `pile-up' at zero problems with Maximum Likelihood Estimation (MLE) of such models. Namely, MLE is likely to generate higher `pile-up' at zero frequencies than MUE when the initial conditions of the state vector are unknown and need to be estimated, and when the true `signal-to-noise ratio' is very small (close to zero). Since holston.etal:2017 do not estimate the initial conditions of the state vector, but instead use very tightly specified prior values, and because their MUEs of the `signal-to-noise ratio' are everything else but very small in the context of MUE, it seems highly unlikely a priori that MLE should generate higher `pile-up' at zero probabilities than MUE. From \cites{stock.watson:1998} simulation results we know that MLE (with a diffuse prior) is substantially more efficient than MUE when the ` signal-to-noise ratio' is not extremely small. MLE should thus be preferred as an estimator.

For reasons of completeness, I provide a comprehensive description of \cites{holston.etal:2017} Stage 1 model and their first stage MUE implementation in \Sref{sec:S1}. As in the Stage 2 model, I show algebraically that their MUE procedure cannot recover the ratio $\sigma _{g}/\sigma _{y^{\ast }}$ from $\lambda _{g}$ because the error term in the first difference of the constructed trend variable $y_{t}^{\ast }$ in the first stage model depends on the real interest rate, as well as `other factor' $z_{t}$ and trend growth $g_{t}$. This means that when the long-run standard deviation from the MUE\ procedure is constructed, it will not only equal $\sigma _{y^{\ast }}$ as required, but also depend on $\sigma _{z}$, $ \sigma _{g}$, as well as the long-run standard deviation of the real rate. Rewriting a simpler version of the Stage 1 model in local level model form also fails to identify the ratio of interest $\sigma _{g}/\sigma _{y^{\ast }} $ from MUE of $\lambda _{g}$. The inability to recover the ratio $\sigma _{g}/\sigma _{y^{\ast }}$ from the first stage model thus appears to be a broader issue highlighting the unsuitability of MUE in this context. This section also illustrates that it is empirically unnecessary to use MUE\ to estimate $\sigma _{g}$ in the first stage model since MLE does not lead to `pile-up' at zero problems with $\sigma _{g}$, neither in the local level model nor in the local linear trend (or unobserved component) model form. Estimating $\sigma _{g}$ directly by MLE\ in the second and third stages confirms this result, yielding in fact larger point estimates than implied by the first stage MUE of $\lambda _{g}$ obtained from \cites{holston.etal:2017} procedure. Readers not interested in the computational intricacies and nuances of the Stage 1 model may skip this section entirely, and only refer back to it as needed for clarification of later results. The key contribution of this paper relates to the correct estimation of $\lambda _{z}$ in \cites{holston.etal:2017} Stage 2 model and its impact on the natural rate of interest through `other factor' $ z_{t}$.

MUE of $\lambda _{z}$ based on the correctly specified Stage 2 model suggests that there is no role for `other factor' $z_{t}$ in this model and given this data.\footnote{ This result is inline with the MLE based estimates of $\sigma _{z}$. Furthermore, these results also carry over to the Euro Area, Canadian and U.K. estimates of $z_{t}$ which are not reported here, but will be made available on the author's webpage.} This brings the focus back to (the estimates of) trend growth in this model. \cites{holston.etal:2017} estimates give the impression that trend growth has markedly slowed since the global financial crisis, particularly in the immediate aftermath of the crisis. In panels (b) and (c) of (ref), I\ show plots of \cites{holston.etal:2017} estimates of $g_{t}$ together with a few simple alternative ones (annualized GDP growth is superimposed in panel (b)). Trend growth is severely underestimated from 2009:Q3 onwards. From the robust (median) estimates of average GDP growth over the various expansion periods shown in (ref), trend growth is only approximately 25 basis points lower at 2.25% since 2009:Q3 than over the pre financial crisis expansion from 2002:Q1 to 2007:Q4.\footnote{ GDP\ growth is close to being serially uncorrelated over the last two expansion periods, with low variances.} Survey based 10 year-ahead expectations of annualized real GDP\ growth plotted in (ref) and (ref) also suggest that trend growth remained stable (these plots are discussed further in \Sref{sec:other}). The key point to take away from this discussion is that \cites{holston.etal:2017} (one sided) Kalman Filter based estimate of $g_{t}$ is excessively `pulled down' by the large decline in GDP during the financial crisis, and this strongly and adversely effects the estimate of trend growth for many periods after the crisis.

The rest of the paper is organised as follows. In \Sref{sec:model}, \cites{holston.etal:2017} structural model of the natural rate of interest is described. \Sref{sec:MUE} gives a concise background to \cites{stock.watson:1998} Median Unbiased Estimation. In \Sref{sec:HLW}, I provide a detailed description of the Stage 1 and Stage 2 models, and report the results of the full Stage 3 model estimates. Some additional issues with the model are discussed in \Sref{sec:other}, and \Sref{sec:conclusion} concludes the study.

\stepcounter{section} \Osubsection*{A.\arabic{section}. {Holston, Laubach and Williams' (2017) Model }} \addcontentsline{toc}{subsection}{\currentname}

\citeallauthors{holston.etal:2017} use the following `structural' model to estimate the natural rate of interest:\footnote{ In what follows, I use the same notation as in holston.etal:2017 (see equations 3 to 9 on pages S61 to S63) to facilitate a direct comparison. Also note that this model builds on an earlier specification of laubach.williams:2003, where trend growth $g_{t}$ is scaled by another parameter $c$, and where also a stationary AR(2) process for the `other factor' $z_{t}$ was considered in addition to the $I(1)$ specification in (\oldref{z}).} \bsq\vsp[-0]

\esq where $y_{t}$ is 100 times the (natural) log of real GDP, $y_{t}^{\ast } $ is the permanent or trend component of GDP, $\tilde{y}_{t}$ is its cyclical component, $\pi _{t}$ is annualized quarter-on-quarter PCE inflation, and $\pi _{t-2,4}=\left( \pi _{t-2}+\pi _{t-3}+\pi _{t-4}\right) /3$. The real interest rate $r_{t}$ is computed as:

where expected inflation is constructed as:

and $i_{t}$ is the exogenously determined nominal interest rate, the federal funds rate.

The natural rate of interest $r_{t}^{\ast }$ is computed as the sum of trend growth $g_{t}$ and `other factor' $z_{t}$, both of which are $I(1)$ processes. The real interest rate gap is defined as $\tilde{r} _{t}=(r_{t}-r_{t}^{\ast })$. The error terms $\varepsilon _{t}^{\ell },\forall \ell =\{\pi ,\tilde{y},y^{\ast }\hsp[-1],g,z\}$ are assumed to be $ i.i.d$ normal distributed, mutually uncorrelated, and with time-invariant variances denoted by $\sigma _{\ell }^{2}$. Notice from (\oldref{AS}) that inflation is restricted to follow an integrated AR(4) process. From the description of the data, we can see that the nominal interest rate $i_{t}$ as well as inflation $\pi _{t}$ are defined in annual or annualized terms, while output, and hence the output gap, trend and trend growth in output are defined at a quarterly rate. Due to this measurement mismatch, holston.etal:2017 adjust the calculation of the natural rate in their code so that trend growth $g_{t}$ is scaled by 4 whenever it enters equations that relate it to annualized variables. The natural rate is thus factually computed as $r_{t}^{\ast }=4g_{t}+z_{t}$.\footnote{ This generates some confusion when working with the model, as it is not clear whether the estimated $z_{t}$ factor is to be interpreted at an annual or quarterly rate.} In the descriptions that follow, I\ will use the annualized $4g_{t}$ trend growth rate whenever it is important to highlight a result or in some of the algebraic derivations, and will leave the equations in (\oldref{eq:hlw}) as in holston.etal:2017 otherwise for ease of comparability.

holston.etal:2017 argue that due to `pile-up' at zero problems with Maximum Likelihood (ML) estimation of the variances of the innovation terms $\varepsilon _{t}^{g}$ and $\varepsilon _{t}^{z}$ in (\oldref{eq:hlw}), estimates of $\sigma _{g}^{2}$ and $\sigma _{z}^{2}$ are \textquotedblleft likely to be biased towards zero \textquotedblright\ (page S64). To avoid such `pile-up' at zero problems, they employ Median Unbiased Estimation (MUE)\ of stock.watson:1998 in two preliminary steps --- Stage 1 and Stage 2 --- to get estimates of what they refer to as `signal-to-noise ratios' defined as $\lambda _{g}=\sigma _{g}/\sigma _{_{y^{\ast }}}$ and $ \lambda _{z}=a_{r}\sigma _{z}/\sigma _{\tilde{y}}$. In Stage 3, the remaining parameters of the full model in (\oldref{eq:hlw}) are estimated, conditional on the median unbiased estimates $\hat{\lambda}_{g}$ and $\hat{ \lambda}_{z}$ obtained in Stages 1 and 2, respectively.

In the above description, I\ intentionally differentiate between the `signal-to-noise ratio' terminology of holston.etal:2017 and the one used in harvey:1989 and in the broader literature on state-space models and exponential smoothing, where the signal-to-noise ratio would be defined as $\sigma _{y^{\ast }}/\sigma _{\tilde{y}}$ or $ \left( \sigma _{g}/\sigma _{\tilde{y}}\right) $ from the relations in (\oldref{eq:hlw}).\footnote{ As noted on page 337 in harvey:2006, the signal-to-noise ratio \textquotedblleft plays the key role in determining how observations should be weighted for prediction and signal extraction.\textquotedblright} To be more explicit, in the context of the classic local level model of muth:1960:\bsq

\esq the signal-to-noise ratio is computed as $\sigma _{\eta }/\sigma _{\varepsilon }$. In the extended version of the model in (\oldref{eq:ll}) known as the local linear trend model: \bsq

\esq two signal-to-noise ratios, namely $\sigma _{\eta }/\sigma _{\varepsilon }$ and $\sigma _{\zeta }/\sigma _{\varepsilon }$, can be formed.\footnote{ The processes $\varepsilon _{t},\eta _{t}$ and $\zeta _{t}$ are uncorrelated white noise. These two state-space formulations are described in more detail in Chapters 2 and 4 of harvey:1989. harvey:1989 also shows how to derive their relation to simple and double exponential smoothing models.} Note here that the model of holston.etal:2017 in (\oldref{eq:hlw}) is essentially an extended and more flexible version of the local linear trend model in (\oldref{eq:llt}). Referring to $\lambda _{g}=\sigma _{g}/\sigma _{_{y^{\ast }}}$ as a signal-to-noise ratio as holston.etal:2017 do is thus rather misleading, since it would correspond to $\sigma _{\zeta }/\sigma _{\eta }$ in the local linear trend model in (\oldref{eq:llt}), which has no relation to the traditional signal-to-noise ratio terminology of harvey:1989 and others in this literature.\footnote{ Readers familiar with the hodrick.prescott:1997 (HP) filter will recognize that the local linear trend model in (\oldref{eq:llt}) --- with the extra `smoothness' restriction $\sigma _{\eta }=0$ --- defines the state-space model representation of the HP\ filter, where the square of the inverse of the signal-to-noise ratio ($\sigma _{\varepsilon }^{2}/\sigma _{\zeta }^{2}$ in (\oldref{eq:llt}) or equivalently $\sigma _{\tilde{y} }^{2}/\sigma _{g}^{2}$ in (\oldref{eq:hlw})) is the HP\ filter smoothing parameter that is frequently set to $1600$ in applications involving quarterly GDP data.}

Before the three stage procedure of holston.etal:2017 is described, I outline in detail how \cites{stock.watson:1998} median unbiased estimator is implemented, what normalization assumptions it imposes, and how look-up tables for the construction of the estimator are computed. I also include a replication of \cites{stock.watson:1998} empirical estimation of trend growth of U.S. real GDP per capita. Although the section that follows below may seem excessively detailed, long, and perhaps unnecessary, the intention here is to provide the reader with an overview of how median unbiased estimation is implemented, what it is intended for, and when one can expect to encounter `pile-up' at zero problems to materialize. Most importantly, it should highlight that `pile-up' at zero is not a problem in the general sense of the word, but rather only a nuisance in situations when it is necessary to distinguish between very small variances and ones that are zero.

\stepcounter{section} \Osubsection*{A.\arabic{section}. {\cites{stock.watson:1998} Median Unbiased Estimation }} \addcontentsline{toc}{subsection}{\currentname}

stock.watson:1998 proposed Median Unbiased Estimation (MUE) in the general setting of Time Varying Parameter (TVP) models. TVP models are commonly specified in a way that allows their parameters to change gradually or smoothly over time. This is achieved by defining the parameters to evolve as driftless random walks (RWs), with the variances of the innovation terms in the RW equations assumed to be small. One issue with Kalman Filter based ML estimation of such models is that estimates of these variances can frequently `pile-up' at zero when the true error variances are `very' small, but nevertheless, non-zero.\footnote{ See the discussion in Section 1 of stock.watson:1998 for additional motivation and explanations. As the title of \cites{stock.watson:1998} paper suggests, MUE was introduced for \textquotedblleft coefficient variance estimation in TVP models\textquotedblright when this variance is expected to be small.}

stock.watson:1998 show simulation evidence of `pile-up' at zero problems with Kalman Filter based ML estimation in Table 1 on page 353 of their paper. In their simulation set-up, they consider the following data generating process for the series $GY_{t}$:\footnote{ See their GAUSS\ files TESTCDF.GSS and ESTLAM.GSS for details on the data generating process, which are available from Mark Watson's homepage at \url{http://www.princeton.edu/ mwatson/ddisk/tvpci.zip}. }\bsq

\esq where $\varepsilon _{t}$ and $\eta _{t}$ are drawn from $i.i.d.$ standard normal distributions, $\beta _{00}$ is initialized at 0, and the sample size is held fixed at $T=500$ observations, using $5000$ replications. The $\lambda $ values that determine the size of the variance of $\Delta \beta _{t}$ are generated over a grid from 0 to 30, with unit increments.\footnote{ To be precise, in their GAUSS\ code, stock.watson:1998 use a range from 0 to 80 for $\lambda $, with finer step sizes for lower $\lambda $ values (see, for instance, the file TESTCDF.GSS). That is, $\lambda $ is a sequence between 0 to 30 with increments of 0.25, then 0.5 unit increments from 30 to 60, and unit increments from 60 to 80. In Tables 1 to 3 of their paper, results are reported for $\lambda $ values from 0 up to $ 30 $ only, with unit increments.} Four median unbiased estimators relying on four different structural break test statistics are compared to two ML estimators. The first ML estimator, referred to as the maximum profile likelihood estimator (MPLE), treats the initial state vector as an unknown parameter to be estimated. The second estimator, the maximum marginal likelihood estimator (MMLE), treats the initial state vector as a Gaussian random variable with a given mean and variance. When the variance of the integrated part of the initial state vector goes to infinity, MMLE produces a likelihood with a diffuse prior.

How one treats the initial condition in the Kalman Filter recursions matters substantially for the `pile-up' at zero problem with MLE. This fact has been known, at least, since the work of shephard.harvey:1990. \footnote{ On page 340, shephard.harvey:1990 write to this: \textquotedblleft \ldots we show that the results for the fixed and known start-up and the diffuse prior are not too different. However, in Section 4 we demonstrate that the sampling distribution of the ML estimator will change dramatically when we specify a fixed but unknown start-up procedure.\textquotedblright Their Tables II and III quantify how much worse the ML estimator that attempts to estimate the initial condition in the local level model performs compared to MLE with a diffuse prior.} The simulation results reported in Table 1 on page 353 in stock.watson:1998 show that `pile-up' at zero frequencies are considerably lower when MMLE with a diffuse prior is used than for MPLE, which estimates the initial state vector. For instance, for the smallest considered non-zero population value of $\lambda =1$, which implies a standard deviation of $\Delta \beta _{t}$ ($\sigma _{\Delta \beta }$ henceforth)\ of $\lambda /T=1/500=0.002$, MMLE produces an at most $14\ $ percentage points higher `pile-up' at zero frequency than MUE (ie., $ 0.60$ or $60\%$ for MMLE versus $0.46$ or $46\%$ for MUE based on the quandt:1960 Likelihood Ratio, henceforth QLR, structural break test statistic).\footnote{ The four different MUEs based on the different structural break tests appear to perform equally well.} For MPLE, this frequency is $45$ percentage points higher at $0.91$ ($91\%$). At $\lambda =5$ ($\sigma _{\Delta \beta }=0.01$) and $\lambda =10$ ($\sigma _{\Delta \beta }=0.02$), these differences in the `pile-up' at zero frequencies reduce to $11$ and $4$ percentage points, respectively, for MMLE, but remain still sizeable for MPLE. At $ \lambda =20$ ($\sigma _{\Delta \beta }=0.04$), the \emph{`pile-up'} at zero problem disappears nearly entirely for MMLE\ and MUE, with \emph{`pile-up'} frequencies dropping to $2$ and $1$ percentage points, respectively, for these two estimators, staying somewhat higher at $7$ percentage points for MPLE.

Using MUE instead of MLE to mitigate `pile-up' at zero problems comes, nevertheless, at a cost; that is, a loss in estimator efficiency whenever $\lambda $ (or $\sigma _{\Delta \beta }$) is not very small. From Table 2 on page 353 in stock.watson:1998, which shows the asymptotic relative efficiency of MUE (and MPLE) relative to MMLE, it is evident that for true $\lambda $ values of $10$ or greater $(\sigma _{\Delta \beta }\geq 0.02)$, the 4 different MUEs yield asymptotic relative efficiencies (AREs) as low as 0.65 (see the results under the $L$ and MW columns in Table 2).\footnote{ The QLR structural break test seems to be the most efficient among the MUEs, yielding the highest AREs across the various MUE implementations.} This means that MMLE only needs $65\%$ of MUE's sample size to achieve the same probability of falling into a given null set. Only for very small values of $ \lambda \leq 4$ $(\sigma _{\Delta \beta }\leq 0.008)$ are the AREs of MUE\ and MMLE of a similar magnitude, ie., close to 1, suggesting that both estimators achieve approximately the same precision.

Three important points are to be taken away from this review of the simulation results reported in stock.watson:1998. First, with MLE, `pile-up' at zero frequencies are substantially smaller when the initial state vector is treated as a known fixed quantity or when a diffuse prior is used, which is the case with MMLE (but not with MPLE). Second, `pile-up' at zero frequencies of MMLE\ are at most $4$ percentage points higher than those of MUE once $\lambda \geq 10$ ($\sigma _{\Delta \beta }=0.02$). Third, MUE can be considerably less efficient than MMLE, in particular for `larger' values of $\lambda \geq 10$ ($\sigma _{\Delta \beta }=0.02$). This suggests that MLE\ with a diffuse prior should be preferred whenever MUE\ based estimates of $\lambda $ (or $\sigma _{\Delta \beta }$) are `large' enough to indicate that `pile-up' at zero problems are unlikely to materialize.

To provide the reader with an illustration of how MUE is implemented, and how its estimates compare to the two maximum likelihood based procedures (MPLE and MMLE), I replicate the empirical example in Section 4 of stock.watson:1998 which provides estimates of trend growth of U.S. real GDP per capita over the period from 1947:Q2 to 1995:Q4. Note that trend growth in GDP is one of the two components that make up the real natural rate $r_{t}^{\ast }$ in holston.etal:2017. It is thus beneficial to illustrate the implementation of MUE in this specific context, rather than in the more general setting of time varying parameter models.

Median unbiased estimation of U.S. trend growth

stock.watson:1998 use the following specification to model the evolution of annualized trend growth in real per capita GDP for the U.S., denoted by $GY_{t}$ below:\footnote{ That is, $GY_{t}=400\Delta \ln (\text{real per capita GDP}_{t})$, where $ \Delta $ is the first difference operator (see Section 4 on page 354 in stock.watson:1998). I again follow their notation as closely as possible for comparability reasons.} \bsq

\esq where $a(L)$ is a (`stationary') lag polynomial with all roots outside the unit circle, $\lambda $ is the parameter of interest, $T$ is the sample size, and $\eta _{t}$ and $\varepsilon _{t}$ are two uncorrelated disturbance terms, with variances $\sigma _{\eta }^{2}$ and $\sigma _{\varepsilon }^{2}$, respectively. The growth rate of per capita GDP is thus composed of a stationary component $u_{t}$ and a random walk component $ \beta _{t}$ for trend growth. stock.watson:1998 set $a(L)$ to a $ 4^{th}$ order lag-polynomial, so that $u_{t}$ follows an AR(4) process. The model in (\oldref{eq:sw98}) can be recognized as the local level model of muth:1960 defined earlier in (\oldref{eq:ll}), albeit with the generalisation that $u_{t}$ follows an AR(4)\ process, rather than white noise. Being in the class of local level models means that the estimate of trend growth will be an exponentially weighted moving average (EWMA) of $GY_{t}$.\footnote{ stock.watson:1998 offer a discussion of the rationale behind the random walk specification of trend growth in $GY_{t}$ in the second paragraph on the left of page 355. Without wanting to get into a technical discussion, one might want to view the random walk specification of trend growth $\beta _{t}$ as a purely statistical tool to allow for a slowly changing mean, rather than interpreting trend growth as an $I(1)$ process.}

It is important to highlight here that \cites{stock.watson:1998} discussion of the theoretical results of the estimator in Sections $2.2-2.3$ of their paper emphasizes that MUE of $\lambda $ in the model in (\oldref{eq:sw98}) is only possible with the \textquotedblleft normalisation $\mathbf{D} =1 $\textquotedblright . They write at the top of page 351\ (right column): \textquotedblleft Henceforth, when $k=1$, we thus set $ \mathbf{D}=1$. When $\mathbf{X}_{t}=1$, under this normalization, $\lambda $\ is $T$\textit{\ times the ratio of the long-run standard deviation of }$\Delta \beta _{t}$\textit{\ to the long run standard deviation of }$u_{t}$\textit{.}\textquotedblright \footnote{ The parameter $k$ here refers to the column dimension of regressor vector $ \mathbf{X}_{t}$. When $k=1$, then only a model with an intercept is fitted, ie., $\mathbf{X}_{t}$ contains only a unit constant and no other regressors.} Denoting the long-run standard deviation of a stochastic process by $\bar{ \sigma}(\cdot )$, this means that

or alternatively, expressed in signal-to-noise ratio form as used by holston.etal:2017:

where $\bar{\sigma}(u_{t})=\sigma _{\varepsilon }/a(1)$ since $u_{t}$ follows a stationary AR(4) process, $a(1)=(1-\sum_{i=1}^{4}a_{i})$, and $ \bar{\sigma}(\Delta \beta _{t})=\sigma _{\Delta \beta }$ due to $\eta _{t}$ being $i.i.d.$, yielding further the relation $\sigma _{\Delta \beta }=(\lambda /T)\sigma _{\eta }$. As a result of the identifying \textquotedblleft normalization $\mathbf{D}=1$\textquotedblright\ of MUE, (\oldref{eq:s2n}) implies that $\sigma _{\eta }=\sigma _{\varepsilon }/a(1)$. That is, the long-run standard deviation of the stationary component $u_{t}$ is equal to the standard deviation of the trend growth innovations $\eta _{t}$.

stock.watson:1998 write on page 354: \textquotedblleft Table 3 is a lookup table that permits computing median unbiased estimates, given a value of the test statistic. The normalization used in Table 3 is that $ \mathbf{D}=1$, and users of this lookup table must be sure to impose this normalization when using the resulting estimator of $\lambda $. \textquotedblright\ Moreover, the numerical results that are reported in\ Section 3, which is appropriately labelled \textquotedblleft Numerical Results for the univariate Local-Level Model\textquotedblright , are obtained from simulations that employ the local level model of (\oldref{eq:tvp_sim}) as the data generating process (see \cites{stock.watson:1998} GAUSS\ programs ESTLAM.GSS, TESTCDF.GSS, and LOOKUP.GSS in the tvpci.zip file that accompanies their paper). These numerical results do not only include the simulations regarding \emph{ `pile-up'} at zero frequencies reported in Table 1, asymptotic power functions plotted in Figure 1, or the AREs provided in Table 2 of stock.watson:1998, but also the look-up tables for the construction of the median unbiased estimator of $\lambda $ in Table 3. It must therefore be kept in mind that these look-up table values are valid only for the univariate local level model, or for models that can be (re-)written in \textit{\textquotedblleft local level form\textquotedblright }.

(ref) below reports the replication results of Tables 4 and 5 in stock.watson:1998.\footnote{ All computations are implemented in Matlab, using their GDP growth data provided in the file DYPC.ASC. Note that I also obtained look-up table values based on a finer grid of $\lambda $ values from their original GAUSS file LOOKUP.GSS (commenting out the lines if (lamdat[i,1] .<= 30) .and (lamdat[i,1]-floor(lamdat[i,1]) .== 0); in LOOKUP.GSS to list look-up values for the entire grid of $\lambda $ 's considered), rather than those listed in Table 3 on page 354 of their paper, where the grid is based on unit increments in $\lambda $ from 0 to 30. I further changed the settings in the tolerance on the gradient in their maximum likelihood (maxlik) library routine to _max_GradTol = 1e-08 and used the printing option format /rd 14,14 for a more precise printing of all results up to 14 decimal points. Lastly, there is a small error in the construction of the lag matrix in the estimation of the AR(4) model in file \texttt{TST_GDP1.GSS} (see lines 40 to 47). The first column in the \texttt{w} matrix is the first lag of the demeaned per capita trend growth series, while columns 2 to 4 are the second to fourth lags of the raw, that is, not demeaned per capita trend growth series. Correcting this leads to mildly higher, yet still insignificant, point estimates of all $ \sigma_{\Delta\beta}$. For instance, the point estimate of $ \sigma_{\Delta\beta}$ based on \cites{nyblom:1989} $L$ statistic yields 0.1501, rather than 0.1303, but remains still statistically insignificant, with the lower value of the confidence interval being 0. To exactly replicate the results in \cites{stock.watson:1998}, I compute the lag matrix as they do.} Columns one and two in the top half of (ref) show test statistics and $p-$values of the four structural break tests that are considered: $i)$ \cites{nyblom:1989} $L$ test, $ii)$ \cites{andrews.ploberger:1994} mean Wald (MW) test, $iii)$ \cites{andrews.ploberger:1994} exponential Wald (EW) test, and $iv)$ \cites{quandt:1960} Likelihood ratio (QLR) test, together with corresponding $p-$values.

As a reminder, the MW, EW and QLR tests are chow:1960 type structural break tests, which test for a structural break in the unconditional mean of a series at a given or known point in time. chow:1960 break tests require a partitioning of the data into two sub-periods. When the break date is unknown, these tests are implemented by rolling through the sample. To be more concrete, denote by $\mathcal{Y}_{t}$ the series to be tested for a structural break in the unconditional mean. Let the dummy variable $ D_{t}(\tau )=1$ if $t>\tau ,$ and $0$ otherwise, where $\tau =\{\tau _{0},\tau _{0}+1,\tau _{0}+2,\ldots ,\tau _{1}\}$ is an index (or sequence) of grid points between endpoints $\tau _{0}$ and $\tau _{1}$. As is common in this literature, stock.watson:1998 set these endpoints at the $ 15^{th}$ and $85^{th}$ percentiles of the sample size $T$, that is, $\tau _{0}=0.15T$ and $\tau _{1}=0.85T$.\footnote{ To be precise, $\tau _{0}$ is computed as $\mathtt{floor(0.15\ast T)}$ and $ \tau _{1}$ as $T-\tau _{0}$. Also, it is standard practice in the structural break literature to trim out some upper/lower percentiles of the search variable to avoid having too few observations at the beginning or at the end of the sample in the 0 and 1 dummy regimes created by $D_{t}(\tau )$. In fact, the large sample approximation of the distribution of the QLR test statistic depends on $\tau _{0}$ and $\tau _{1}$. stock.watson:2011 write to this on page 558: \textquotedblleft For the large-sample approximation to the distribution of the QLR statistic to be a good one, the sub-sample endpoints, $\tau _{0}$ and $\tau _{1}$, cannot be too close to the beginning or the end of the sample.\textquotedblright\ Employing endpoints other than the $15^{th}$ upper/lower percentile values used by stock.watson:1998 in the simulation of the look-up table for $ \lambda $ is thus likely to affect the values provided in Table 3 of stock.watson:1998, due to the endpoints' influence on the distribution of the structural break test statistics.} For each $\tau \in \lbrack \tau _{0},\tau _{1}]$, the following regression of $\mathcal{Y}_{t}$ on an intercept term and $D_{t}(\tau )$ is estimated:

and the $F$ statistic (the square of the $t-$statistic) on the point estimate $\hat{\zeta}_{1}$ is constructed. The sequence $\{F(\tau )\}_{\tau =\tau _{0}}^{\tau _{1}}$ of $F$ statistics is then utilized to compute the MW, EW and QLR structural break test statistics needed in the implementation of MUE. These are calculated as:\bsq

\esq where $N_{\tau }$ denotes the number of grid points in $\tau $. \cites{nyblom:1989} $L$ test statistic is computed without sequentially partitioning the data via the sum of squared cumulative sums of $\mathcal{Y} _{t}$. More specifically, let $\hat{\mu}_{\mathcal{Y}}$ denote the sample mean of $\mathcal{Y}_{t}$, $\hat{\sigma}_{\mathcal{Y}}^{2}$ the sample variance of $\mathcal{Y}_{t}$, and $\mathcal{\tilde{Y}}_{t}=\mathcal{Y}_{t}- \hat{\mu}_{\mathcal{Y}}$ the demeaned $\mathcal{Y}_{t}$ process. \cites{nyblom:1989} $L$ statistic is then constructed as:

where $\vartheta _{t}$ is the scaled cumulative sum of $\mathcal{\tilde{Y}} _{t}$, ie., $\vartheta _{t}=T^{-1/2}\sum_{s=1}^{t}\mathcal{\tilde{Y}}_{s} $.

Median unbiased estimates of $\lambda $ based on \cites{stock.watson:1998} look-up tables are reported in column 3 of (ref), followed by respective 90% confidence intervals (CIs) in square brackets. The last two columns show estimates of $\sigma _{\Delta \beta }$ computed as $\hat{ \sigma}_{\Delta \beta }=\hat{\lambda}/T\times \hat{\sigma}_{\varepsilon }/ \hat{a}(1)$, with 90% CIs also in square brackets.\ In the bottom half of (ref), MLE\ and MUE based parameter estimates of the model in (\oldref{eq:sw98}) are reported. The columns under the MPLE and MMLE headings show, respectively, MLE based results when the initial state vector is estimated and when a diffuse prior is used. The diffuse prior for the $I(1)$ element of the state vector is centered at 0 with a variance of $10^{6}$. The next two columns under the headings MUE(0.13) and MUE(0.62) report parameter estimates of the model in (\oldref{eq:sw98}) with $\sigma _{\Delta \beta }$ held fixed at its MUE$\ $point estimate of $0.13$ and upper 90% CI value of $0.62$, respectively. The last column under the heading SW.GAUSS lists the corresponding MUE(0.13) estimates obtained from running \cites{stock.watson:1998} GAUSS\ code as reference values.\footnote{ See the results reported in Table 5 on page 354 in stock.watson:1998, where nevertheless only two decimal points are reported. MPLE and MMLE are also replicated accurately to 6 decimal points.}

As can be seen from the results in (ref), consistent with the `pile-up' at zero problem documented in the simulations in stock.watson:1998 (and also shephard.harvey:1990), the MPLE\ estimate of $\sigma _{\Delta \beta }$ goes numerically to zero (up to 11 decimal points), while MMLE produces a `sizeable' point estimate of $\sigma _{\Delta \beta }$ of $0.044$. Although stock.watson:1998 (and also\ I) do not report a standard error for $\hat{\sigma}_{\Delta \beta }$ in the tables containing the estimation results, the estimate of $\mathrm{ stderr}(\hat{\sigma}_{\Delta \beta })$ is $0.1520$, suggesting that $\hat{ \sigma}_{\Delta \beta }$ is very imprecisely estimated.\footnote{stock.watson:1998 compute standard errors for the remaining MMLE\ parameters (see column three in the upper part of Table 5 on page 354 in their paper. They write in the notes to Table 5: \textquotedblleft Because of the nonnormal distribution of the MLE of $\lambda $, the standard error for $\sigma _{\Delta \beta }$ is not reported .\textquotedblright\ Evidently, `testing' the null hypothesis of $ \sigma _{\Delta \beta }=0$ using a standard $t-$ratio does not make any sense statistically. Nevertheless, $\hat{\sigma}_{\Delta \beta }$ is very imprecisely estimated, and highly likely to be \emph{`very'} close to zero. The MMLE log-likelihood function with the restriction $\sigma _{\Delta \beta }=0$ is $-547.5781$, while the (unrestricted) MMLE is $-547.4805$, with the difference between the two being very small of about $0.10$.} From the MUE\ results reported in the first column of the top half of (ref) it is evident that all 4 structural break tests yield confidence intervals for $\lambda $ and hence also $\sigma _{\Delta \beta }$ that include zero. Thus, even when using MUE as the \textit{`preferred'} estimator, one would conclude that $\hat{\lambda}$ and $\hat{\sigma}_{\Delta \beta }$ are \textit{ not} statistically different from zero.

An evident practical problem with the use of \cites{stock.watson:1998} MUE\ is that the 4 different structural break tests can produce vastly different point estimates of $\lambda $. This is clearly visible from (ref), where the 4 tests yield $\lambda $ estimates with an implied $ \hat{\sigma}_{\Delta \beta }$ range between $0.0250$ (for QLR)$\ $and $ 0.1303 $ (for $L$). From the simulation results in stock.watson:1998 we know that all 4 tests seem to behave equally well in the `pile-up' at zero frequency simulations (see Table 1 in stock.watson:1998). However, the QLR test performed `best' in the efficiency results, producing the largest (closes to 1)\ asymptotic relative efficiencies in Table 2 of stock.watson:1998. Analysing these results in the context of the empirical estimation of trend growth, the most accurate MUE estimator based on the QLR\ structural break test produces an estimate of $\sigma _{\Delta \beta }$ that is 5 times smaller than the largest one based on the $L$ structural break test, with the MMLE\ estimate of $\sigma _{\Delta \beta } $ being approximately double the size of the QLR estimate.

To provide a visual feel of how different the MLE\ and MUE based estimates of U.S. trend growth are, I show plots of the smoothed estimates in (ref) (these correspond to Figures 4 and 3 in stock.watson:1998). The top panel displays the MPLE, MMLE, MUE(0.13), and MUE(0.62) estimates together with a 90%\ CI of the MMLE estimate (shaded area), as well as a dashed yellow line that shows \cites{stock.watson:1998} GAUSS code based MUE(0.13) estimate for reference.\ The plot in the bottom panel of (ref) superimposes the actual $GY_{t}$ series to portray the variability in the trend growth estimates relative to the variation in the data from which these were extracted.\footnote{ Notice from the top panel of (ref) that there are four different estimates of $\lambda $, and thus four $\hat{\sigma}_{\Delta \beta }$. Rather then showing smoothed trend estimates for all four of these, I follow stock.watson:1998 and only show estimates based on \cites{nyblom:1989} $L$ statistic, which has the largest $\lambda ~$ estimate, and hence also $\sigma _{\Delta \beta }$.} The $y-$axis range is set as in Figures 4 and 3 in stock.watson:1998. As can be seen from (ref), there is only little variability in the MLE based trend growth estimates, with somewhat more variation from MUE(0.13). Nonetheless, all three trend growth point estimates stay within the 90% error bands of MMLE. Moreover, the plots in (ref) confirm the lack of precision of MUE. Trend growth could be anywhere between a constant value of about $1.8\%$ ($\hat{\beta}_{00}$ from MPLE), which is a flat line graphically when $\sigma _{\Delta \beta }$ is held fixed at its lower $90\%$\ CI value of 0, and a rather volatile series which produces a range between nearly $4.5\%$ in 1950 and less than $0.5\%$ in 1980 when $ \sigma _{\Delta \beta }$ is set at its upper $90\%$\ CI value of 0.62.

Given the previous results and discussion, one could argue that the statistical evidence in support of any important time varying trend growth in real U.S.\ GDP per capita is rather weak in this model and data set. As a robustness check and in the context of a broader replication of the time varying trend growth estimates of stock.watson:1998, I obtain real GDP per capita data from the Federal Reserve Economic Data (FRED2) database and re-estimate model. These results are reported in (ref) , which is arranged in the same way as (ref) (only the last column with heading SW.GAUSS is removed). The sample period is again from 1947:Q2 to 1995:Q4, using an AR(4)\ model to approximate $u_{t}$ in (\oldref{eq:sw1}).\footnote{ The results using an ARMA($2,2$) model for $u_{t}$ instead are qualitatively the same.} From (ref) it is clear that not only do the two MLE based estimates of $\sigma _{\Delta \beta }$ yield point estimates that are numerically equal to zero, but so do all 4 MUEs. Hence, trend growth may well be constant. More importantly, it demonstrates that MUE\ can also lead to zero estimates of $\sigma _{\Delta \beta }$ and that there is nothing unusual about that.\footnote{ I show later that the Stage 2 MUE procedure of holston.etal:2017 is incorrectly implemented and based on a misspecified Stage 2 model. Once this is corrected, the Stage 2 $\lambda _{z}$ that one obtains is very close to zero, resulting in the full model MLE\ and MUE\ estimates being very similar. }

Before I proceed to describe how the three stage procedure of holston.etal:2017 is implemented, a brief procedural description of \cites{stock.watson:1998} MUE\ that lists the main steps needed to replicate the results reported in (ref) and (ref) is provided below.

\stepcounter{section} \Osubsection*{A.\arabic{section}. {Three stage estimation procedure of holston.etal:2017 }} \addcontentsline{toc}{subsection}{\currentname}

holston.etal:2017 employ MUE in two preliminary stages that are based on restricted versions of the full model in (\oldref{eq:hlw}) to obtain estimates of the `signal-to-noise ratios' $\lambda _{g}=\sigma _{g}/\sigma _{_{y^{\ast }}}$ and $\lambda _{z}=a_{r}\sigma _{z}/\sigma _{\tilde{y}}$. These ratios are then held fixed in Stage 3 of their procedure, which produces estimates of the remaining parameters of the model in (\oldref{eq:hlw}). In order to conserve space in the main text, I provide all algebraic details needed for the replication of the three individual stages in the \hyperref[appendix]{Appendix}, which includes also some additional discussion as well as R-Code extracts to show the exact computations. In the results that are reported in this section, I have used their R-Code from the file \href{https://www.newyorkfed.org/medialibrary/media/research/economists/williams/data/HLW_Code.zip} {HLW_Code.zip} made available on Willams' website at the New York Fed to (numerically) accurately reproduce their results.\footnote{ Williams' website at the Federal Reserve Bank of New York is at: \url{https://www.newyorkfed.org/research/economists/williams/pub}. Their R-Code is available from the website: \url{https://www.newyorkfed.org/medialibrary/media/research/economists/williams/data/HLW_Code.zip} . The weblink to the file with their real time estimates is: \url{https://www.newyorkfed.org/medialibrary/media/research/economists/williams/data/Holston_Laubach_Williams_real_time_estimates.xlsx} . Note here that all my results exactly match their estimates provided in the \href{https://www.newyorkfed.org/medialibrary/media/research/economists/williams/data/Holston_Laubach_Williams_real_time_estimates.xlsx} {Holston_Laubach_Williams_real_time_estimates.xlsx} file in Sheet 2017Q1.} The sample period that I cover ends in 2017:Q1.\ The beginning of the sample is the same as in holston.etal:2017. That is, it starts in 1960:Q1, where the first 4 quarters are used for initialisation of the state vector, while the estimation period starts in 1961:Q1.

holston.etal:2017 adopt the general state-space model (SSM) notation of hamilton:1994 in their three stage procedure. The SSM is formulated as follows:\footnote{ The state-space form that they use is described on pages 9 to 11 of their online appendix that is included with the R-Code HLW_Code.zip file from Williams' website at the New York Fed. Note that I use exactly the same state-space notation to facilitate the comparison to holston.etal:2017 , with the only exception being that I include one extra selection matrix term $\mathbf{S}$ in front of $\boldsymbol{\epsilon }_{t}$ in (\oldref{eq:RQ}) as is common in the literature to match the dimension of the state vector to $ \boldsymbol{\epsilon }_{t}$ when there are identities due to lagged values. I also prefer not to transpose the system matrices $\mathbf{A}$\ and \ $\mathbf{H}$ in (\oldref{eq:RQ}), as it is not necessary and does not improve the readability.}

where we can define $\boldsymbol{\epsilon }_{t}=\mathbf{S}\boldsymbol{ \varepsilon }_{t}$, so that $\mathrm{Var}(\boldsymbol{\epsilon }_{t})= \mathrm{Var}(\mathbf{S}\boldsymbol{\varepsilon }_{t})=\mathbf{SWS}^{\prime }= \mathbf{Q}$ to make it consistent with the notation used in holston.etal:2017. The (observed) measurement vector is denoted by $ \mathbf{y}_{t}$ in (\oldref{eq:RQ}), $\mathbf{x}_{t}$ is a vector of exogenous variables, $\mathbf{A}$, $\mathbf{H}$ and $\mathbf{F}$ are conformable system matrices, $\boldsymbol{\xi }_{t}$ is the latent state vector, $\mathbf{S}$ is a selection matrix, and the notation $\mathsf{MNorm} \left( \boldsymbol{\mu },\boldsymbol{\Sigma }\right) $ denotes a multivariate normal random variable with mean vector $\boldsymbol{\mu }$ and covariance matrix $\boldsymbol{\Sigma }$. The disturbance terms $\boldsymbol{ \nu }_{t}$ and $\boldsymbol{\varepsilon }_{t}$ are serially uncorrelated, and the (individual) covariance matrices $\mathbf{R}$ and $\mathbf{W}$ are assumed to be diagonal matrices, implying zero correlation between the elements of the measurement and state vector disturbance terms. The measurement vector $\mathbf{y}_{t}$ in (\oldref{eq:RQ}) is the same for all three stages and is defined as $\mathbf{y}_{t}=[y_{t},~\pi _{t}]^{\prime }$, where $y_{t}$ and $\pi _{t}$ are the log of real GDP\ and annualized PCE inflation, respectively, as defined in \Sref{sec:model}. The exact form of the remaining components of the SSM\ in (\oldref{eq:RQ}) changes with the estimation stage that is considered, and is described in detail either in the text below or in the \hyperref[appendix]{Appendix}.

As I have emphasized in the description of MUE in \Sref{sec:MUE}, the simulation results of stock.watson:1998 show that `pile-up' at zero frequencies for MLE are not only a function of the size of the variance of $\Delta \beta _{t}=(\lambda /T)\eta _{t}$ (or alternatively $\lambda $), but also depend critically on whether the initial condition of the state vector is estimated or not. Now holston.etal:2017 do not estimate the initial condition of the state vector in any of the three stages that are implemented. Instead, they apply the HP filter to log GDP data with the smoothing parameter set to $36000$ to get a preliminary estimate of $y_{t}^{\ast }$ and trend growth $g_{t}$ (computed as the first difference of the HP\ filter estimate of $y_{t}^{\ast }$) using data from 1960:Q1 onwards. `Other factor' $z_{t}$ is initialized at 0. \footnote{ See the listing in \hyperref[R:stage3]{R-Code\nref{R:stage3} {1}} in the \hyperref[Rcode]{A.6 R-Code Snippets} section of the \hyperref[appendix]{Appendix}, which shows the first 122 lines of their R-file rstar.stage3.R. Line 30 shows the construction of the initial state vector as $\boldsymbol{\xi }_{00}=[y_{0}^{\ast },y_{-1}^{\ast },y_{-2}^{\ast },g_{-1},g_{-2},z_{-1},z_{-2}]^{\prime }$ where subscripts $ [0,-1,-2]$ refer to the time periods 1960:Q4, 1960:Q3, and 1960:Q2, respectively. In terms of their R-Code, we have: xi.00 <- c(100* g.pot[3:1],100*g.pot.diff[2:1],0,0), where \texttt{g.pot} is the HP filtered trend and \texttt{g.pot.diff} is its first difference, ie., trend growth, with the two zeros at the end being the initialisation of $ z_{t}$. This yields the following numerical values: [806.45, 805.29, 804.12, 1.1604, 1.1603, 0, 0]. The same strategy is also used in the first two stages (see their R-files \texttt{rstar.stage1.R} and \texttt{rstar.stage2.R} ).} This means that $\boldsymbol{\xi }_{00}$ has known and fixed quantities in all three stages.\ Given the simulation evidence provided in Table 1 on page 353 in stock.watson:1998, one may thus expect a priori \emph{`pile-up'} at zero frequencies of MLE (without estimation of the initial conditions) to be only marginally larger than those of MUE, especially for everything but very small values of $\lambda $.

Also, holston.etal:2017 determined the covariance matrix of the initial state vector in an unorthodox way. Even though every element of the state vector $\boldsymbol{\xi }_{t}$ in all three estimation stages is an $ I(1)$ variable, they do not employ a diffuse prior on the state vector. Instead, the covariance matrix is determined with a call to the function calculate.covariance.R (see the code snippet in \hyperref[R:covar]{R-Code\nref{R:covar} {2}} for details on this function, and also lines 66, 84, and 88, respectively, in their R-files rstar.stage1.R, rstar.stage2.R, and rstar.stage3.R, with line 88 in rstar.stage3.R also shown on the second page of the code snippet in \hyperref[R:stage3]{R-Code\nref{R:stage3} {1}}). To summarize what this function does, consider the Stage 1 model, which is estimated with a call to rstar.stage1.R. The function \texttt{ calculate.covariance.R }first sets the initial covariance matrix to $0.2$ times a three dimensional identity matrix $\mathbf{I}_{3}$. Their procedure then continues by using data from 1961:Q1 to the end of the sample to get an estimate of $\sigma _{y^{\ast }}^{2}$ from the Stage 1 model. Lastly, the initial covariance matrix $\mathbf{P}_{00}$ to be used in the \textit{`} \emph{final}\textit{'} estimation of the Stage 1 model is then computed as: \bsq

\esq with $\mathbf{\hat{Q}}$ a $(3\times 3)$ dimensional zero matrix with element $(1,1)$ set to $\hat{\sigma}_{y^{\ast }}^{2}=0.27113455739$ from the initial run of the Stage 1 model. What this procedure effectively does is to set $\mathbf{P}_{00}$ to the first time period's predicted state covariance matrix, given an initial state covariance matrix of $0.2\times \mathbf{I} _{3} $ and the estimate $\hat{\sigma}_{y^{\ast }}^{2}$, where $\hat{\sigma} _{y^{\ast }}^{2}$ was obtained by MLE\ and the Kalman Filter using $ 0.2\times \mathbf{I}_{3}$ as the initial state covariance. This way of initialising $\mathbf{P}_{00}$ is rather circular, as it fundamentally presets $\mathbf{P}_{00}$ at $0.2\times \mathbf{I}_{3}$.\footnote{ In footnote 6 on page S64 in holston.etal:2017 (and also in the description of the calculate.covariance.R file), they write: \textquotedblleft We compute the covariance matrix of these states from the gradients of the likelihood function.\textquotedblright\ Given the contents of the R-Code, it is unclear how and if this was implemented.}

When the state vector contains $I(1)$ variables, it is not only standard practice to use a diffuse prior, but it is highly recommended. For instance, harvey:1989 writes to this on the bottom of page 121: \textquotedblleft When the transition equation is non-stationary, the unconditional distribution of the state vector is not defined. Unless genuine prior information is available, therefore, the initial distribution of $\boldsymbol{\alpha }_{0}$\ must be specified in terms of a diffuse or non-informative prior.\textquotedblright\ (emphasis added, $ \boldsymbol{\alpha }_{t}$ is the state-vector in Harvey's notation). It is not clear why holston.etal:2017 do not use a diffuse prior.\footnote{ In an earlier paper using a similar model for the NAIRU, laubach:2001 discusses the use of diffuse priors. laubach:2001 writes on page 222:\ \textquotedblleft The most commonly used approach in the presence of a nonstationary state is to integrate the initial value out of the likelihood by specifying an (approximately) diffuse prior.\textquotedblright\ He then proceeds to describe an alternative procedure that can be implemented by using: \textquotedblleft a few initial observations to estimate the initial state by GLS, and use the covariance matrix of the estimator as initial value for the conditional covariance matrix of the state.\textquotedblright\ The discussion is then closed with the statement: \textquotedblleft This is the first approach considered here. Because this estimate of the initial state and its covariance matrix are functions of the model parameters, under certain parameter choices the covariance matrix may be ill conditioned. The routines then choose the diffuse prior described above as default.\textquotedblright \ Thus even here, the diffuse prior is the "\emph{safe}" default option. Note that their current procedure does not use: \textquotedblleft \emph{a few initial observations to estimate the initial state}\textquotedblright , but the same sample of data that are used in the final model, ie., with data beginning in 1961:Q1.} However, one may conjecture that it could be due to their preference for reporting Kalman Filtered (one-sided) rather than the more efficient Kalman Smoothed (two-sided) estimates of the latent state vector $\boldsymbol{\xi }_{t}$ which includes trend growth $g_{t}$ and \emph{ `other factor'} $z_{t}$ needed to construct $r_{t}^{\ast }$.\footnote{ Note that Filtered estimates of $g_{t}$, $z_{t}$ and thus also $r_{t}^{\ast } $ are very volatile at the beginning of the sample period (until about 1970) when $\mathbf{P}_{00}$ is initialized with a diffuse prior.}

As a final point in relation to the probability of `pile-up' at zero problems arising due to small variances of the state innovations, and hence the rationale for employing MUE rather than MLE\ in the first place, one can observe from the size of the $\sigma _{g}$ and $\sigma _{z}$ estimates for the U.S. reported in Table 1 on page S60 in holston.etal:2017 that these are rather `large' at $0.122$ and $0.150$, respectively. The simulation results in Table 1 in stock.watson:1998 show that `pile-up' at zero frequencies drop to $0.01$ for both, MMLE\ and MUE, when the true population value of $\lambda $ is $30$ ($\sigma _{\Delta \beta }=0.06$). Given the fact that holston.etal:2017 do not estimate the initial value of the state vector, and that their median unbiased estimates are about two times larger than $0.06$, it seems highly implausible that `pile-up' at zero problems should materialize with a higher probability for MLE than for MUE.

Stage 1 model

\cites{holston.etal:2017} first stage model takes the following restricted form of the full model presented in equation (\oldref{eq:hlw}):\footnote{ See \hyperref[sec:AS1]{Section A.1} in the \hyperref[appendix]{Appendix} for the exact matrix expressions and expansions of the first stage SSM. Note that one key difference of \cites{holston.etal:2017} SSM specification described in equations (\oldref{AS1:m}) and (\oldref{AS1:s}) in the \hyperref[appendix] {Appendix} is that the expansion of the system matrices for the Stage 1 model does not include the drift term $g$ in the trend specification in (\oldref{S1d}), so that $y_{t}^{\ast }$ follows a random walk without drift. Evidently, such a specification cannot match the upward trend in the GDP data. To resolve this mismatch, holston.etal:2017 `detrend' output $y_{t}$ in the estimation (see \hyperref[sec:AS1]{Section A.1} in the \hyperref[appendix]{Appendix} which describes how this is done and also shows snippets of their R-Code).}\bsq

\esq where the vector of Stage 1 parameters to be estimated is:

To be able to distinguish the disturbance terms of the full model in (\oldref{eq:hlw}) from the ones in the restricted Stage 1 model in (\oldref{eq:stag1}) above, I have placed a ring $(\mathring{\phantom{y}})$ symbol on the error terms in (\oldref{S1c}) and (\oldref{S1d}). These two disturbance terms from the restricted model are defined as:

and

From the relations in (\oldref{S1eps_ystar0}) and (\oldref{S1eps_ytilde0}) it is clear that, due to the restrictions in the Stage 1 model, the error terms $ \mathring{\varepsilon}_{t}^{\tilde{y}}$ and $\mathring{\varepsilon} _{t}^{y^{\ast }}$ in (\oldref{eq:stag1}) will not be uncorrelated anymore, since $ \mathrm{Cov}(\mathring{\varepsilon}_{t}^{\tilde{y}},\mathring{\varepsilon} _{t}^{y^{\ast }})=-\tfrac{a_{r}}{2}4\sigma _{g}^{2}$ given the assumptions of the full model in (\oldref{eq:hlw}).\ The separation of trend and cycle shocks in this formulation of the Stage 1 model is thus more intricate, as both shocks will respond to one common factor, the missing $g_{t-1}$.

In the implementation of the Stage 1 model, holston.etal:2017 make two important modelling choices that have a substantial impact on the $ \boldsymbol{\theta }_{1}$ parameter estimates, and thus also the estimate of the `signal-to-noise ratio' $\lambda _{g}$ used in the later stages. The first is the tight specification of the prior variance of the initial state vector $\mathbf{P}_{00}$ discussed in the introduction of this section. The second is a lower bound restriction on $b_{y}$ in the inflation equation in (\oldref{S1b}) ($b_{y}\geq 0.025$ in the estimation). The effect of these two choices on the estimates of the Stage 1 model parameters are shown in (ref) below. The left block of the estimates in (ref) (under the heading `HLW Prior') reports four sets of results where the state vector was initialized using their values for $\boldsymbol{ \xi }_{00}$ and $\mathbf{P}_{00}$. The first column of this block (HLW.R-File) reports estimates from running \cites{holston.etal:2017} R-Code for the first stage model. These are reported as reference values. The second column ($b_{y}\geq 0.025$) shows my replication of \cites{holston.etal:2017} results using the same initial values for parameter vector $\boldsymbol{\theta }_{1}$ in the optimisation routine and also the same lower bound constraint on $b_{y}$. The third column (Alt.Init.Vals) displays the results I obtain when a different initial value for $b_{y}$ is used, with the lower bound restriction $b_{y}\geq 0.025$ still in place. The fourth column ($b_{y} $ Free)\ reports results when the lower bound constraint on $b_{y}$ is removed.\footnote{ To find the initial values for $\boldsymbol{\theta }_{1}$, holston.etal:2017 apply the HP filter to GDP\ to obtain an initial estimate of the cycle and trend components of GDP. These estimates are then used to find initial values for (some of) the components of parameter vector $\boldsymbol{\theta }_{1}$ by running OLS\ regressions of the HP cycle estimate on two of its own lags (an AR(2) essentially), and by running regressions of inflation on its own lags and one lag of the HP cycle. Interestingly, although readily available, rather than taking the coefficient on the lagged value of the HP cycle in the initialization of $ b_{y}$, which yields a value of $0.0921$, holston.etal:2017 use the lower bound value of $0.025$ for $b_{y}$ as the initial value. In the optimisation, this has the effect that the estimate for $b_{y}$ is effectively stuck at $0.025$, although it is not the global optimum in the restricted model, which is at $b_{y}=0.097185$ (see also the values of the log-likelihood function reported in the last row of (ref)).} The right block in (ref) shows parameter estimates when a diffuse prior for $\boldsymbol{\xi }_{t}$ is used, where $\mathbf{P}_{00}$ is set to $10^{6}$ times a three dimensional identity matrix, with the left and right columns showing, respectively, the estimates with and without the lower bound restriction on $b_{y}$ imposed.

Notice initially from the first two columns in the left block of (ref) that their numerical results are accurately replicated up to 6 decimal points. From these results we also see that the lower bound restriction on $b_{y}$ is binding. holston.etal:2017 set the initial value for $b_{y}$ at $0.025$, and there is no movement away from this value in the numerical routine. Specifying an alternative initial value for $b_{y}$ , which is determined in the same way as for the remaining parameters in $ \boldsymbol{\theta }_{1}$, leads to markedly different estimates, while removing the lower bound restriction on $b_{y}$ all together results in the ML estimate of $b_{y}$ to converge to zero. Evidently, these three scenarios yield also noticeably different values for ${\hat{\sigma}}_{y^{\ast }}$, that is, values between $0.4190$ and $0.6177$. The diffuse prior based results (with and without the lower bound restriction) in the right block of (ref) show somewhat less variability in ${\hat{\sigma}} _{y^{\ast }}$, but affect the persistence of the cycle variable $\tilde{y} _{t}$ in the model, with the smallest AR(2) lag polynomial root being 1.1190 when $b_{y}\geq 0.025$ is imposed, while it is only 1.0251 and thus closer to the unit circle when $b_{y}$ is left unrestricted.

As a final comment, there is only little variation in the likelihoods of the different estimates that are reported in the respective left and right blocks of (ref). For instance, the largest difference in log-likelihoods is obtained from the diffuse prior results shown in the right block of (ref). If we treat the lower bound as a restriction, a Likelihood Ratio (LR) test of the null hypothesis of the difference in these likelihoods being zero yields $ -2(-536.9803-(-535.9596))=2.0414$, which, with one degree of freedom has a $ p-$value of $0.1531$ and cannot be rejected at conventional significance levels. Hence, there is only limited information in the data to compute a precise estimate of $b_{y}$. This empirical fact is known in the literature as a `flat Phillips curve'.\footnote{ That the output gap is nearly uninformative for inflation (forecasting) once structural break information is conditioned upon --- regardless of what measure of the output gap is used or whether it is combined as an ensemble from multiple measures --- is shown in buncic.muller:2017 for the U.S. and for Switzerland.}

Given the Stage 1 estimate $\skew{0}\boldsymbol{\hat{\theta}}_{1}$, holston.etal:2017 use the following steps to implement median unbiased estimation of their `signal-to-noise ratio' $\lambda _{g}=\sigma _{g}/\sigma _{y^{\ast }}$.

(ref) shows the range of $\lambda _{g}$ estimates computed from the five sets of $\skew{0}\boldsymbol{\hat{\theta}}_{1}$ values reported in (ref), using all four structural break tests of stock.watson:1998. (ref) is arranged in the same format as (ref), again showing \cites{holston.etal:2017} estimates of $\lambda _{g}$ obtained from running their R-Code in the first column of the left block for reference. As can be seen from (ref), the range of $\hat{\lambda}_{g}$ values one obtains from \cites{holston.etal:2017} MUE procedure is between $ 0 $ to $0.08945$ (if only the three structural break tests implemented by holston.etal:2017 are considered, and up to $0.09419$ if the $L$ statistic is computed as well. Note that this range is not due to statistical uncertainty, but simply due to the choice of structural break test, which prior for $\mathbf{P}_{00}$ is used, and whether the lower bound constraint on $b_{y}$ is imposed. Since these estimates determine the relative variation in trend growth through the magnitude of $\sigma _{y^{\ast }}$, they have a direct impact not only on the variation in the permanent component of GDP, but also on the natural rate of interest through the ratio $\lambda _{g}=\sigma _{g}/\sigma _{y^{\ast }}$ utilized in the later stages of the three step procedure of holston.etal:2017.

\cites{holston.etal:2017} rational for MUE in Stage 1

Comparing the MUE procedure that holston.etal:2017 implement to the one by stock.watson:1998, it is evident that they are fundamentally different. Instead of rewriting the true model of interest in local level form to make it compatible with \cites{stock.watson:1998} look-up tables, holston.etal:2017 instead formulate a restricted Stage 1 model that not only sets $a_{r}$ in the output gap equation to zero, but also makes the awkward assumption that trend growth is constant when computing the `preliminary' estimate of $y_{t}^{\ast }$.

The rationale behind \cites{holston.etal:2017} implementation of MUE is as follows. Suppose we observe trend $y_{t}^{\ast }$. Then, a local level model for $\Delta y_{t}^{\ast }$ can be formulated as:\bsq

\esq where $\Delta y_{t}^{\ast }$, $g_{t}$ and $\varepsilon _{t}^{y^{\ast }}$ are the analogues to $GY_{t},\beta _{t}$ and $u_{t}$, respectively, in \cites{stock.watson:1998} MUE in (\oldref{eq:sw98}), with $\varepsilon _{t}^{y^{\ast }}$ in (\oldref{hlw_ll1a}), nonetheless, assumed to be $i.i.d.$ rather than an autocorrelated AR(4) process as $u_{t}$ in (\oldref{eq:sw1}). Under \cites{stock.watson:1998} assumptions, MUE of the local level model in (\oldref{eq:hlwLL1}) yields $\lambda _{g}=\lambda /T$ defined as:

where $\bar{\sigma}(\cdot )$ denotes again the long-run standard deviation, and the last equality in (\oldref{hlw_rationale1}) follows due to $\varepsilon _{t}^{y^{\ast }}$ and $\varepsilon _{t}^{g}$ assumed to be uncorrelated white noise processes.

Since $\Delta y_{t}^{\ast }$ is not observed, holston.etal:2017 replace it with the Kalman Smoother based estimate $\Delta \hat{y} _{t|T}^{\ast }$ obtained from the restricted Stage 1 model in (\oldref{eq:stag1}) . To illustrate what impact this has on their MUE procedure, let $ a_{y}(L)=(1-a_{y,1}L-a_{y,2}L^{2})$ and $a_{r}(L)=\tfrac{a_{r}}{2}(L+L^{2})$ denote two lag polynomials that capture the dynamics in the output gap $ \tilde{y}_{t}$ and the real rate cycle $\tilde{r}_{t}=(r_{t}-r_{t}^{\ast })=(r_{t}-4g_{t}-z_{t})$, respectively. Also, define $\psi (L)=a_{y}(L)^{-1}a_{r}(L)$ and $\psi (1)=a_{r}/(1-a_{y,1}-a_{y,2})$. The output gap equation of the true (full) model in (\oldref{eq:hlw}) can then be written compactly as:

Observed output, and trend and cycle are related by the identity

This relation, together with (\oldref{hlw_ll1a}) and (\oldref{s2dyb}), can be written as:

\qquad

Because the data $\Delta y_{t}$ are fixed, any restriction imposed on the $ \Delta \tilde{y}_{t}$ process translates directly into a misspecification of the right hand side of (\oldref{eq:mis1}); the $\Delta y_{t}^{\ast }$ term. In the Stage 1 model, $a_{r}$ is restricted to zero. For the relation in (\oldref{eq:mis1}) to balance, $\Delta y_{t}^{\ast }$ effectively becomes:\footnote{ Note that we need to formulate a local level model for trend growth as in (\oldref{eq:hlwLL1}) to be able to apply the MUE\ framework of stock.watson:1998. To arrive at (\oldref{hlw_ll2a}), add $ [a_{y}(L)^{-1}a_{r}(L)\Delta \tilde{r}_{t}]$ to both sides of (\oldref{eq:mis1}). The ring $(\mathring{\phantom{y}})$ symbol on $\mathring{\nu}_{t}^{y\ast }$ highlights again that it is obtained from the restricted model. } \bsq

\esq where

\cites{holston.etal:2017} implementation of MUE relies on the (constructed) local level model relations from the restricted Stage 1 model in (\oldref{S1LLfalse}) and requires us to evaluate the ratio of the long-run standard deviations of $\varepsilon _{t}^{g}$ and $\mathring{\nu}_{t}^{y\ast }$:

Evidently, $\varepsilon _{t}^{g}$ in (\oldref{hlw_ll2b}) has not changed, so the numerator of the `signal-to-noise ratio' in (\oldref{s2n_s1}) is still $ \bar{\sigma}(\varepsilon _{t}^{g})=\sigma _{g}\ $, due to $\varepsilon _{t}^{g}$ being an $i.i.d.$ process. However, the term $\mathring{\nu} _{t}^{y\ast }$ in (\oldref{hlw_ll2a}) is not uncorrelated white noise anymore. Moreover, the long-run standard deviation $\bar{\sigma}(\mathring{\nu} _{t}^{y\ast })$ in the denominator of (\oldref{s2n_s1}) now also depends on the (long-run) standard deviation of $\psi (L)\Delta \tilde{r}_{t}$, and will be equal to $\sigma _{y^{\ast }}$ if and only if $a_{r}=0$ in the empirical data.\footnote{ If monetary policy is believed to be effective in cyclical aggregate demand management, then $a_{r}$ cannot be 0 and one would not have formulated the main model of interest assuming that $a_{r}$ is different from zero (viz, negative). Also, this restriction cannot be enforced in the data.}

To see what the long-run standard deviation of $\mathring{\nu}_{t}^{y\ast }$ looks like, assume for simplicity that $\varepsilon _{t}^{y^{\ast }}$ and $ \Delta \tilde{r}_{t}$ are uncorrelated, so that the long-run standard deviation calculation of $\mathring{\nu}_{t}^{y\ast }$ can be broken up into a part involving $\varepsilon _{t}^{y^{\ast }}$ and another part involving $ \psi (L)\Delta \tilde{r}_{t}$, where the latter decomposes as:

Assuming that the shocks $\{\varepsilon _{t}^{g},\varepsilon _{t}^{z}\}$ are uncorrelated with the (change in the) real rate $\Delta r_{t}$, the long-run standard deviation of $\psi (L)\Delta \tilde{r}_{t}$ can be evaluated as:

since $\varepsilon _{t}^{g}$ and $\varepsilon _{t}^{z}$ are uncorrelated in the model. Because the nominal rate $i_{t}$ is exogenous, it will not be possible to say more about the first term on the right hand side of (\oldref{LRv}) unless we assume some time series process for $\Delta r_{t}$. Suppose that $ r_{t}$ follows a random walk, so that $\Delta r_{t}=\varepsilon _{t}^{r}$, with $\mathrm{Var}(\varepsilon _{t}^{r})=\sigma _{r}^{2}$. Then $\bar{\sigma} \left( \psi (L)\Delta \tilde{r}_{t}\right) =a_{r}/(1-a_{y,1}-a_{y,2})\left[ \sigma _{r}+4\sigma _{g}+\sigma _{z}\right] $, and we obtain $\bar{\sigma}( \mathring{\nu}_{t}^{y\ast })=\sigma _{y^{\ast }}+a_{r}/(1-a_{y,1}-a_{y,2}) \left[ \sigma _{r}+4\sigma _{g}+\sigma _{z}\right] $. The MUE ratio in (\oldref{s2n_s1}) based on the restricted Stage 1 model yields:

Thus, \cites{holston.etal:2017} implementation of MUE in Stage 1 cannot recover the `signal-to-noise ratio' of interest $\frac{ \sigma _{g}}{\sigma _{y^{\ast }}}$ from $\lambda _{g}$.

Note here that the autocorrelation pattern in $\mathring{\nu}_{t}^{y\ast }$ is also reflected in the $\Delta \hat{y}_{t|T}^{\ast }$ series which is used as the observable counterpart to $\Delta y_{t}^{\ast }$ in (\oldref{hlw_ll2a}). That is, $\Delta \hat{y}_{t|T}^{\ast }$ has a significant and sizeable AR(1) coefficient of $-0.2320$ (standard error $\approx 0.0649$). Inline with Step $(i)$ of \cites{stock.watson:1998} implementation of MUE (the GLS\ step), one would thus need to AR(1)\ filter the constructed $\Delta \hat{y} _{t|T}^{\ast }$ series used in the local level model before implementing the structural break tests. Accounting for this autocorrelation patter in $\Delta \hat{y}_{t|T}^{\ast }$ leads to very different $\lambda _{g}$ point estimates (see (ref), which is arranged in the same way as the top half of (ref), with the last column showing $\lambda _{g}=\lambda /T$ rather than $\sigma _{g}$ to be able to compare these to column one of (ref)).

Rewriting the Stage 1 model in local level model form

One nuisance with the Stage 1 model formulation of holston.etal:2017 in (\oldref{eq:stag1}) is that trend growth is initially assumed to be constant to compute a first estimate of $y_{t}^{\ast }$. This estimate is then used to construct the empirical counterpart of $\Delta y_{t}^{\ast }$ to which MUE is applied.

A more coherent way to implement MUE in the context of the Stage 1 model is to rewrite the local linear trend model in local level form. To see how this could be done, we can simplify the Stage 1 model by excluding the inflation equation (\oldref{S1b}) and replacing the constant trend growth equation in (\oldref{S1d}) with the original trend and trend growth equations in (\oldref{y*}) and (\oldref{g}). Since the specification of the full model in (\oldref{eq:hlw}) assumes that the error terms $\varepsilon _{t}^{\ell },\forall \ell =\{\pi ,\tilde{y} ,y^{\ast }\hsp[-1],g,z\}$ are $i.i.d.$ Normal and mutually uncorrelated, and $\hat{b}_{y}\approx 0$ in the unrestricted Stage 1 model (see the results under the heading `$b_{y}$ Free' in (ref)), this simplification is unlikely to induce any additional misspecification into the model.

The modified Stage 1 model we can work with thus takes the following form: \bsq

\esq where $\mathring{\varepsilon}_{t}^{\tilde{y}}=a_{r}(L)\tilde{r} _{t}+\varepsilon _{t}^{\tilde{y}}$ again due to the restriction of the output gap equation of the full model in (\oldref{eq:hlw}).\footnote{ If the disturbance term $\mathring{\varepsilon}_{t}^{\tilde{y}}$ is $i.i.d.$ , then the model in (\oldref{Stage1:mod}) can be recognized as \cites{clark:1987} Unobserved Component (UC) model. However, $\mathring{\varepsilon}_{t}^{ \tilde{y}}$ is not $i.i.d.$ and instead follows a general ARMA\ process with non-zero autocovariances, which are functions of $\sigma _{g}^{2}$, $\sigma _{z}^{2}$, the autocovariances of inflation $\pi _{t}$, as well as the exogenously specified interest rate $i_{t}$. To see this, recall from \Sref{sec:model} that the real interest rate gap $\tilde{r}_{t}$ is defied as $\tilde{r}_{t}=\left[ i_{t}-\delta (L)\pi _{t}-4g_{t}-z_{t}\right] $, where expected inflation $\pi _{t}^{e}=\delta (L)\pi _{t}$ and $\delta (L)= \tfrac{1}{4}\left( 1+L+L^{2}+L^{3}\right) $, so that we can re-express $ \mathring{\varepsilon}_{t}^{\tilde{y}}$ as:

The product of the two lag polynomials $a_{r}(L)\delta (L)$ in (\oldref{eps_tilde}) yields a $5^{th}$ order lag polynomial for inflation. If $i_{t}$ and $\pi _{t}$ were uncorrelated white noise processes (which they are clearly not), then we would obtain an MA(5)\ process for $\mathring{ \varepsilon}_{t}^{\tilde{y}}$ when $a_{r}$ is non-zero. Since $\pi _{t}$ is modelled as an integrated AR(4), the implied process for $\mathring{ \varepsilon}_{t}^{\tilde{y}}$ is a higher order ARMA\ process, the exact order of which depends on the assumptions one places on the exogenously specified interest rate $i_{t}$. To determine this process exactly is of no material interest here. However, the important point to take away from this is that $\mathring{\varepsilon}_{t}^{ \tilde{y}}$ is autocorrelated and follows a higher order ARMA\ process. Moreover, if $i_{t},\pi _{t},g_{t}$ and $z_{t}$ do not co-integrate, then $\mathring{\varepsilon}_{t}^{\tilde{y} } $ will be an $I(1)$ process.} The local linear trend model in (\oldref{Stage1:mod}) can now be rewritten in local level model form by differencing (\oldref{S1M1a}) and (\oldref{S1M1b}), and bringing $y_{t-1}^{\ast }$ to the left side of (\oldref{S1M1c}) to give the relations:\bsq

\esq

Substituting (\oldref{S1M2b}) and (\oldref{S1M2c}) into (\oldref{S1M2a}) yields the local level model:

where $u_{t}$ is defined as:

with $b(L)\varepsilon _{t}=a_{y}(L)\varepsilon _{t}^{y^{\ast }}+\Delta \mathring{\varepsilon}_{t}^{\tilde{y}}$ on the right hand side of (\oldref{eq:ut_arma}) denoting a general MA process. The $u_{t}$ term in (\oldref{eq:ut_arma}) thus follows a higher order ARMA model. If $a_{r}=0$, then $ \mathring{\varepsilon}_{t}^{\tilde{y}}=\varepsilon _{t}^{\tilde{y}}$ in (\oldref{eps_tilde}) and $\Delta \mathring{\varepsilon}_{t}^{\tilde{y}}=\Delta \varepsilon _{t}^{\tilde{y}}$, which is an integrated MA(1)\ process, so that the right hand side would be the sum of an MA(2)\ and an MA(1), yielding an overall MA(2) for $b(L)\varepsilon _{t}$. With $a_{y}(L)$ being an AR(2) lag polynomial for the cycle component, we would then get an ARMA$ (2,2)$ for $u_{t}$ in (\oldref{eq:ut_arma}). If $a_{r}\neq 0$, then $\Delta \mathring{\varepsilon}_{t}^{\tilde{y}}$ follows a higher order ARMA process. In the empirical implementation of MUE, I follow stock.watson:1998, and use an AR(4) as an approximating model for $u_{t}$.\footnote{ They also considered an ARMA($2,3$) model (see page 355 in their paper). It is well known that higher order ARMA\ models can be difficult to estimate numerically due to potential root cancellations in the AR and MA lag polynomials. Inspection of the autocorrelation and partial autocorrelation functions of $\Delta y_{t}$ indicate that an AR(4) model is more than adequate to capture the time series dynamics of $\Delta y_{t}$. I have also estimated an ARMA$(2,2)$\ model for $\Delta y_{t}\,$, with the overall qualitative conclusions being the same and the quantitative results very similar.}

The relations in (\oldref{S1Mod2b}) to (\oldref{eq:ut_arma}) are now in local level model form to which MUE\ can be applied to as outlined in equations (\oldref{eq:sw98}) to (\oldref{eq:s2n}) in \Sref{subsec:MUE}.\footnote{ I am grateful to James Stock for his email correspondence on this point.} To examine if we can recover the `signal-to-noise ratio' of interest $ \sigma _{g}/\sigma _{y^{\ast }}$ from this MUE procedure, we need to evaluate

In the numerator of (\oldref{eq:S1Lambda}), the term $\bar{\sigma}(\varepsilon _{t}^{g})=\sigma _{g}$ as before. Nevertheless, the denominator term $\bar{ \sigma}(u_{t})=\bar{\sigma}(\varepsilon _{t}^{y^{\ast }}+a_{y}(L)^{-1}\Delta \mathring{\varepsilon}_{t}^{\tilde{y}})\neq \sigma _{y^{\ast }}$. With $ \mathring{\varepsilon}_{t}^{\tilde{y}}=a_{r}(L)\tilde{r}_{t}+\varepsilon _{t}^{\tilde{y}}$, we have:\

where the middle part in (\oldref{LRv2}) (ie., $\psi (L)\Delta \tilde{r}_{t}$) will again be as before in (\oldref{LRv}) and therefore depend on $\Delta r_{t}$, $\varepsilon _{t}^{g}$ and $\varepsilon _{t}^{z}$. Notice here also that even if we knew $a_{r}=0$, so that the middle part in (\oldref{LRv2}) is $0$, there is no mechanism to enforce a zero correlation between $\varepsilon _{t}^{y^{\ast }}$ and $\varepsilon _{t}^{\tilde{y}}$ in the data, because $ u_{t}$ appears in reduced form in the local level model. We would thus need the empirical correlation between $\varepsilon _{t}^{y^{\ast }}$ and $ \varepsilon _{t}^{\tilde{y}}$ to be zero for the long-run standard deviation $\bar{\sigma}(u_{t})$ to equal $\sigma _{y^{\ast }}$ even when the true $ a_{r}=0$. Estimates from the existing business cycle literature suggest that trend and cycle shocks are negatively correlated (see for instance Table 3 in morley.etal:2003, who estimate this correlation to be $-0.9062$, or Table 1 in the more recent study by grant.chan:2017a whose estimate is $-0.87$). I obtain an estimate of $-0.9426$ (see (ref) below).

For completeness, parameter estimates of MUE applied to the local level transformed Stage 1 model defined in (\oldref{Stage1:mod}) are reported in (ref). This table is arranged in the same way as (ref), with all computations performed in exactly the same way as before. The MUE results in the last two columns of the bottom part of the table are based on the exponential Wald (EW) structural break test as used in holston.etal:2017. Overall, these estimates are very similar to \cites{stock.watson:1998} estimates, despite different time periods and GDP data being used. The $\lambda $ (and also $\sigma _{g}$) estimates are not statistically different from 0, and the MMLE $\hat{\sigma}_{g}$ of $0.1062$ is rather sizeable and quite close to the one implied by MUE.

Estimating the local linear trend version of the Stage 1 model

So far, `pile-up' at zero problems were examined in the local level model form which is compatible with MUE. As a last exercise, I estimate the modified Stage 1 model in (\oldref{Stage1:mod}) in local linear trend model form. Two different specifications of the model are estimated. The first assumes all error terms to be uncorrelated. This version is referred to as \cites{clark:1987} UC0 model. The second allows for a non-zero correlation between $\varepsilon _{t}^{y^{\ast }}$ and $\mathring{\varepsilon}_{t}^{ \tilde{y}}$. This version is labelled \cites{clark:1987} UC model. The aim here is to not only examine empirically how valid the zero correlation assumption is and to quantify its magnitude, but also to investigate whether `pile-up' at zero problems materialize more generally in UC\ models. In (ref), the parameter estimates of the two UC models are reported, together with standard errors of the parameter estimates (these are listed under the columns with the heading Std.error).

As can be seen from the estimates in (ref), there exists no evidence of `pile-up' at zero problems with MLE in either of these two UC models.\footnote{ I use a diffuse prior on the initial state vector in the estimation of both UC models, and do not estimate the initial value. This is analogous to MMLE in stock.watson:1998. The input data are $100$ times the log of real GDP.} The estimates of $\sigma _{g}$ from the two UC models are $0.0463$ and $0.0322$, respectively, and are based on quarterly data. Expressed at an annualized rate, they amount to approximately $0.1852$ and $0.1288$, and hence are similar in magnitude to the corresponding MUE based estimates obtained from the transformed model in (ref). Notice also that the correlation between $\mathring{\varepsilon}_{t}^{\tilde{y}}$ and $ \varepsilon _{t}^{y^{\ast }}$ (denoted by $\mathrm{Corr}(\mathring{ \varepsilon}_{t}^{\tilde{y}},\varepsilon _{t}^{y^{\ast }})$ in (ref)) is estimated to be $-0.9426$ ($t-$statistic is approximately $-10$). The magnitude of the $\hat{\sigma}_{y^{\ast }}$ and $ \hat{\sigma}_{\tilde{y}}$ coefficients nearly doubles when an allowance for a non-zero correlation between $\mathring{\varepsilon}_{t}^{\tilde{y}}$ and $ \varepsilon _{t}^{y^{\ast }}$ is made.\footnote{ As is common with UC models, the improvement in the log-likelihood due to the addition of the extra correlation parameter is rather small. Although it is important to empirically capture the correlation between $\mathring{ \varepsilon}_{t}^{\tilde{y}}$ and $\varepsilon _{t}^{y^{\ast }}$ as it affects the trend growth estimate (see (ref)), the overall level of information contained in the data appears to be limited and therefore makes it difficult to decisively distinguish one model over the other statistically. Also, one other aspect of the empirical GDP data that both models fail to capture is the global financial crisis. The level of GDP dropped substantially and in an unprecedented manner. Simply `smoothing' the data to extract a trend as the UC models implicity do may thus not adequately capture this drop in the level of the series.}

(ref) shows plots of the various trend growth estimates from the modified Stage 1 models reported in (ref) and (ref). The plots are presented in the same way as in (ref) earlier, with the (annualized) trend growth estimates from the two UC models superimposed. Analogous to the results in stock.watson:1998, the variation in the MUE based estimates is once again large. Trend growth can be a flat line when the lower 90% CI\ of MUE is considered or rather variable when the upper CI bound is used. Interestingly, the MMLE, Clark UC model (with non-zero $\mathrm{Corr}( \mathring{\varepsilon}_{t}^{\tilde{y}},\varepsilon _{t}^{y^{\ast }})$) and MUE$(\hat{\lambda}_{EW})$ trend growth estimates are very similar visually. More importantly, the effect of restricting $\mathrm{Corr}(\mathring{ \varepsilon}_{t}^{\tilde{y}},\varepsilon _{t}^{y^{\ast }})$ to zero on the trend growth estimate can be directly seen in (ref). The UC0 model produces a noticeably more variable trend growth estimate than the UC model.

Two conclusions can be drawn from this section. Firstly, \cites{holston.etal:2017} implementation of MUE in Stage 1 and the resulting $\lambda _{g}$ estimate cannot recover the `signal-to-noise ratio' of interest $\sigma _{g}/\sigma _{y^{\ast }}$. Secondly, there is no evidence of `pile-up' at zero problems materializing when estimating $\sigma _{g}$ directly by MLE. Replacing $\sigma _{g}$ in $ \mathbf{Q}$ by $\hat{\lambda}_{g}\sigma _{y^{\ast }}$ in the Stage 2 and full model log-likelihood functions (see (\oldref{S2Q}) and (\oldref{Q3})) where $\hat{ \lambda}_{g}$ was obtained from MUE applied to the Stage 1 model is not only unsound but empirically entirely unnecessary.

Stage 2 Model

The second stage model of holston.etal:2017 consists of the following system of equations, which are again a restricted version of the full model in (\oldref{eq:hlw}):\bsq

\esq Given the estimate of $\lambda _{g}$ from Stage 1, the vector of Stage 2 parameters to be estimated by MLE is:\footnote{ See \hyperref[sec:AS2]{Section A.2} in the \hyperref[appendix]{Appendix} for the exact matrix expressions and expansions of the SSM of Stage 2. In the $ \mathbf{Q}$ matrix, $\sigma _{g}$ is replaced by $\hat{\lambda}_{g}\sigma _{y^{\ast }}$, where $\hat{\lambda}_{g}$ is the estimate from the first stage model (see (\oldref{S2Q})). The state vector $\boldsymbol{\xi }_{t}$ is initialized using the same procedure as outlined in (\oldref{eq:P00S1a}) and \fnref{fn:1}, with the numerical values of $\boldsymbol{\xi }_{00}$ and $ \mathbf{P}_{00}$ given in (\oldref{AS2:xi00}) and (\oldref{AS2:P00}).}

As in the first stage model in (\oldref{eq:stag1}), I again use the ring symbol $( \mathring{\phantom{y}})$ on the disturbance terms in (\oldref{S2:ytilde}) and (\oldref{S2:ystar}) to distinguish them from the $i.i.d.$ error terms of the full model in (\oldref{eq:hlw}).

Examining the formulation of the Stage 2 model in (\oldref{eq:stag2}) and comparing it to the full model in (\oldref{eq:hlw}), it is evident that holston.etal:2017 make two `misspecification' choices that are important to highlight. First, they include $g_{t-2}$ instead of $g_{t-1}$ in the trend equation in (\oldref{S2:ystar}), so that the $\mathring{\varepsilon} _{t}^{y^{\ast }}$ error term is in fact:\footnote{holston.etal:2017 only report the $\mathbf{\mathbf{Q}}$ matrix in their documentation, which is a diagonal matrix and takes the form given in (\oldref{S2Q}). In \hyperref[sec:AS2] {Section A.2} of the \hyperref[appendix]{Appendix}, I\ show how this matrix is obtained. In \hyperref[sec:AS21]{Section A.2.1}, the correct Stage 2 model state-space form is provided, applying the same `trick' as used in the Stage 3 state-space model specification. The two $\mathbf{\mathbf{Q}}$ matrices are listed in (\oldref{S2Q}) and (\oldref{S2Qcorrect}).}

As a result of this, $\mathring{\varepsilon}_{t}^{y^{\ast }}$ in (\oldref{e_ystar}) follows an MA(1)\ process, instead of white noise as $\varepsilon _{t}^{y^{\ast }}$ in (\oldref{y*}). Moreover, due to the $\varepsilon _{t-1}^{g}$ term in (\oldref{e_ystar}), the covariance between the two error terms in (\oldref{S2:ystar}) and (\oldref{S2:g}) is no longer zero, but rather $\sigma _{g}^{2}$. Thus, treating $\mathbf{W}$ in (\oldref{eq:RQ}) as a diagonal variance-covariance matrix in the estimation of the second stage model is incorrect.

Second, holston.etal:2017 do not only add an (unnecessary) intercept term $a_{0}$ to the output gap equation in (\oldref{S2:ytilde}), but they also account for only one lag in trend growth $g_{t}$, and further fail to impose the $a_{g}=-4a_{r}$ restriction in the estimation of $a_{g}$. Due to this, the error term $\mathring{\varepsilon}_{t}^{\tilde{y}}$ in (\oldref{S2:ytilde}) can be seen to consist of the following two components:

where the `desired terms' on the right-hand side of (\oldref{S2:eps_ytilde}) are needed for \cites{holston.etal:2017} implementation of MUE in the second stage, whose logic I will explain momentarily, while the `unnecessary terms' are purely due to the ad hoc addition of an intercept term, changing lag structure on $g_{t}$ and failure to impose the $a_{g}=-4a_{r}$ restriction.

To be consistent with the full model specification in (\oldref{eq:hlw}), the relations in (\oldref{S2:ytilde}) and (\oldref{S2:ystar}) should have been formulated as: \bsq

\esq so that only the two missing lags of $z_{t}$ from (\oldref{S2a:ytilde}) appear in the error term $\mathring{\varepsilon}_{t}^{\tilde{y}}$, specifically:

Such a specification could have been easily obtained from the full Stage 3 state-space model form described in \hyperref[sec:AS3]{Section A.3} in the \hyperref[appendix]{Appendix}, by simply removing the last two row entries of the state vector $\boldsymbol{\xi }_{t}$ in (\oldref{AS3:xi}), and adjusting the $\mathbf{H}$, $\mathbf{F}$, and $\mathbf{S}$ matrices in the state and measurement equations to be conformable with this state vector. This is illustrated in \hyperref[sec:AS21]{Section A.2.1} in the \hyperref[appendix]{ Appendix}. The `correctly specified' Stage 2 model should thus have been:\bsq

\esq

To see why this matters, let us examine how one would implement MUE\ in the Stage 2 model, following again \cites{holston.etal:2017} logic as applied in Stage 1. That is, one would first need to define a local level model involving $z_{t}$ to be in the same format as in (\oldref{eq:sw98}). If we assume for the moment that the true state variables $\tilde{y}_{t}$ and $g_{t}$, as well as parameters $a_{y,1},$ $a_{y,2}$ and $a_{r}$ are known, and we ignore the econometric issues that arise when these are replaced by estimates, then the following local level model from the `correctly specified' Stage 2 model in (\oldref{S2_ytilde0}) can be formed:\bsq

\esq where $a_{y}(L)\tilde{y}_{t}-a_{r}(L)[r_{t}-4g_{t}]$ and $ -a_{r}(L)z_{t} $ in (\oldref{MUE2a}) are the analogues to $GY_{t}$ and $\beta _{t} $ in (\oldref{eq:sw1}), $\varepsilon _{t}^{\tilde{y}}$ corresponds to $u_{t}$ (but is $i.i.d.$ from the full model assumptions in (\oldref{eq:hlw}) rather than an autocorrelated time series process as $u_{t}$ in (\oldref{eq:sw1})), and $ -a_{r}(L)\Delta z_{t}$ and $-a_{r}(L)\varepsilon _{t}^{z}$ are the counterparts to $\Delta \beta _{t}$ and $(\lambda /T)\eta _{t}$ in the state equation in (\oldref{eq:swRW}).\footnote{ To arrive at (\oldref{MUE2b}), simply multiply (\oldref{z}) in the full model by $ -a_{r}(L)$.}

The equations in (\oldref{MUE2}) are now in local level model form suitable for MUE. The Stage 2 MUE procedure implemented on this constructed $ GY_{t}=a_{y}(L)\tilde{y}_{t}-a_{r}(L)[r_{t}-4g_{t}]$ series produces the $ \lambda _{z}=\lambda /T$ ratio corresponding to (\oldref{eq:s2n}), that is: \footnote{ To make this clear, MUE returns an estimate of $\lambda $ by using the look-up table on page 354 in stock.watson:1998 to find the closest matching value of one of the four structural break test statistics defined in (\oldref{eq:breakTests}) and (\oldref{eqL}) which test for a structural break in the unconditional mean of the constructed $GY_{t}$ series by running a dummy variable regression of the form defined in (\oldref{Zt}).}

The last two steps in (\oldref{mue2_ratio}) follow due to $a_{r}(1)=\frac{a_{r}}{2 }(1+1^{2})=a_{r}$ and $\bar{\sigma}(\varepsilon _{t}^{\tilde{y}})=\sigma _{ \tilde{y}}$, with $\bar{\sigma}(\cdot )$ denoting again the long-run standard deviation. The final term in (\oldref{mue2_ratio}) gives \cites{holston.etal:2017} ratio $\lambda _{z}=a_{r}\sigma _{z}/\sigma _{ \tilde{y}}$.\footnote{ In laubach.williams:2003, $\lambda _{z}$ is curiously defined as the ratio $a_{r}\sigma _{z}/(\sigma _{\tilde{y}}\sqrt{2})$ (see page 1064, second paragraph on the right). It is not clear where the extra $\sqrt{2}$ term comes from.} This is the logic behind \cites{holston.etal:2017} implementation of MUE in Stage 2.

However, because holston.etal:2017 define the Stage 2 model in `misspecified' form in (\oldref{eq:stag2}), $\mathring{\varepsilon}_{t}^{ \tilde{y}}$ is no longer simply equal to $-a_{r}(L)z_{t}+\varepsilon _{t}^{ \tilde{y}}$ as needed for the right-hand side of (\oldref{MUE2a}), but now also includes the `unnecessary terms' $\left[ a_{0}+a_{g}g_{t-1}+a_{r}(L)4g_{t}\right] $ (see the decomposition in (\oldref{S2:eps_ytilde})). What effect this has on the Stage 2 MUE\ procedure can be seen by first rewriting $a_{g}g_{t-1}$ as:

where $a_{g}(L)=\frac{a_{g}}{2}(L+L^{2})$. The additional `unnecessary terms' on the right-hand side of (\oldref{S2:eps_ytilde}) become:\vsp[-2]

In \cites{holston.etal:2017} Stage 2 model in (\oldref{eq:stag2}), the constructed local level model takes then the form:\bsq

\esq where $\mathring{\nu}_{t}^{\tilde{y}}$ in (\oldref{S2wrong_a}) is the misspecified analogue to $\varepsilon _{t}^{\tilde{y}}$ in (\oldref{MUE2a}) and is defined as:

As can be seen, the error term $\mathring{\nu}_{t}^{\tilde{y}}$ in (\oldref{S2_nu_ring}) will not be white noise. Moreover, forming the MUE\ $\lambda /T$ ratio from the model in (\oldref{S2wrong}) in the same way as in (\oldref{mue2_ratio}) leads to:

and now requires the evaluation of the long-run standard deviation of $ \mathring{\nu}_{t}^{\tilde{y}}$ in the denominator, which will not be equal to $\sigma _{\tilde{y}}$ as from the `correctly' specified Stage 2 model defined in (\oldref{S2full0}). Note here that, even in the unlikely scenario that $(a_{g}+4a_{r})=0$ in the data, the long-run standard deviation of $\mathring{\nu}_{t}^{\tilde{y}}$ will also depend on $\tfrac{ a_{g}}{2}\sigma _{g}$ because of the $\tfrac{a_{g}}{2}\varepsilon _{t-1}^{g}$ term in $\mathring{\nu}_{t}^{\tilde{y}}$, so that one obtains:

Thus, MUE applied to \cites{holston.etal:2017} `misspecified' Stage 2 model as defined in (\oldref{eq:stag2}) cannot recover the ratio of interest $ \lambda _{z}=a_{r}\sigma _{z}/\sigma _{\tilde{y}}$.\footnote{ If $(a_{g}+4a_{r})\neq 0$, additional $\sigma _{g}$ terms enter the long-run standard deviation in the denominator of $\lambda _{z}$.}

Before I\ discuss in the next section what effect\ the `misspecification' of the Stage 2 model has on \cites{holston.etal:2017} median unbiased estimates of $\lambda _{z}$, I\ report the estimates of the two different Stage 2 models in (ref). The first and second columns show replicated results which are based on \cites{holston.etal:2017} R-Code as well as my own implementation and serve as reference values. In the third column under the heading `MLE$(\sigma _{g})$', $\sigma _{g}$ is estimated directly by MLE together with the other parameters of the model without using $\hat{\lambda}_{g}$ from Stage 1.\footnote{ I use the same initial values for the parameter and the state vector (mean and variance) as in the exact replication of holston.etal:2017. Using a diffuse prior instead leads to only minor differences in the numerical values. The implied $\lambda _{g}$ and $\sigma _{g}$ estimates are shown in brackets and were computed from the `signal-to-noise ratio' relation $ \lambda _{g}=\sigma _{g}/\sigma _{y^{\ast }}.$} The last column under the heading `MLE$(\sigma _{g}).\mathcal{M}_{0}$' reports estimates obtained from the `correctly specified' Stage 2 model defined in (\oldref{S2full0}), where $\sigma _{g}$ is once again estimated directly by MLE.

The results in (ref) can be summarized as follows. First, there exists no evidence of `pile-up' at zero problems materializing when estimating $\sigma _{g}$ directly by MLE; not in the `misspecified' Stage 2 model, nor in the `correctly specified' \ one. This finding is consistent with the earlier results from the first stage. The Stage 2 MLE of $\sigma _{g}$ is in fact nearly $50\%$ larger than the estimate implied by $\hat{\lambda}_{g}$ from MUE in Stage 1. MUE in Stage 1 thus seems to be redundant. Second, the estimate of $ a_{g}$ is about eight times the magnitude of $-a_{r}$, so that $ (a_{g}+4a_{r})\approx 0.3132\neq 0$. Therefore, the ratio in (\oldref{Lz00}) will have additional $\sigma _{g}$ terms in the denominator, making the evaluation of this quantity more intricate. And third, despite the different Stage 2 model specifications, the resulting parameter estimates as well as the log-likelihood values across the three different models in columns two to four of (ref) are very similar. This suggests that, overall, the data are uninformative about the model parameters.\footnote{ These findings also hold when using data for the Euro Area, the U.K., and Canada, but are not reported here.}

Note here that, although the results in (ref) indicate that `misspecifying' the Stage 2 model does not have an important impact on the parameter estimates that are obtained, I show below that it substantially and spuriously amplifies the size of the $\lambda _{z}$ estimate.

\cites{holston.etal:2017} implementation of MUE in Stage 2

Recall again conceptually how MUE in Stage 2 would need to be implemented following the same logic as in Stage 1 before.\ First, one needs to construct an observable counterpart to $GY_{t}$ as given in (\oldref{MUE2a}) from the Stage 2 model estimates. Then, the four structural break tests described in \Sref{subsec:MUE} are applied to test for a break in the unconditional mean of (the AR filtered) $GY_{t}$ series. This corresponds to Step $(ii)$ in \cites{stock.watson:1998} procedural description. Constructing a local level model of the form described in (\oldref{MUE2}) enables us to implement MUE\ to yield the ratio $\lambda /T=\bar{\sigma}(\Delta \beta _{t})/\bar{\sigma} (\varepsilon _{t}^{\tilde{y}})$ as defined in (\oldref{mue2_ratio}).

\cites{holston.etal:2017} implementation of MUE in Stage 2, nonetheless, departs from this description in two important ways. First, instead of using the `correctly specified' Stage 2 model defined in (\oldref{S2full0}), they work with the `misspecified' model given in (\oldref{eq:stag2}). Second, rather than leaving\ the $a_{y,1},$ $a_{y,2},$ $a_{r},$ $a_{g}$ and $ a_{0}$ parameters fixed at their Stage 2 estimates and constructing the observable counterpart to $GY_{t}$ in (\oldref{S2wrong_a}) only once outside the dummy variable regression loop, holston.etal:2017 essentially `re-estimate' these parameters by including the vector $\boldsymbol{ \mathcal{X}}_{t}$ defined in (\oldref{XX}) below as a regressor in the structural break regression in (\oldref{eqS2regs}). For the `misspecified' Stage 2 model, this has the effect of substantially increasing the size and variability of not only the dummy variable coefficients $\hat{\zeta}_{1}$ in (\oldref{eqS2regs}), but also the corresponding $F$ statistics used in the computation of the MW, EW, and \textrm{QLR} structural break tests needed for MUE$\ $of $\lambda _{z}$.

To illustrate how holston.etal:2017 implement MUE in the second stage, I list below the main steps that they follow to compute $\lambda _{z}$ .

In the top and bottom panels of (ref) I show plots of the sequences of $F$ statistics $\{F(\tau )\}_{\tau =\tau _{0}}^{\tau _{1}}$ computed from \cites{holston.etal:2017} `misspecified' Stage 2 model and the `correctly specified' Stage 2 model defined in (\oldref{S2full0}), respectively. Two sets of sequences are drawn in each panel. \footnote{ The same sequence computed from an updated data series up to 2019:Q2 is shown in (ref) in the \hyperref[appendix]{Appendix}.} The first sequence, which I refer to as `time varying $\boldsymbol{ \phi }$' (drawn as a red line in (ref)) is constructed by following \cites{holston.etal:2017} implementation outlined in Steps (I\hsp[.3]) to (\emph{III}\hsp[.3]) above. I call this the \emph{`time varying} $\boldsymbol{\phi }$\textit{'} sequence because the $a_{y,1},$ $a_{y,2},$ $a_{r},$ $a_{g}$ and $a_{0}$ parameters needed to \textit{`construct'} the observable counterpart to $ GY_{t}$ in (\oldref{S2wrong_a}) are effectively \textit{`re-estimated'} for each $ \tau \in \lbrack \tau _{0},\tau _{1}]$ in the dummy variable regression loop due to the inclusion of the extra $\boldsymbol{\mathcal{X}}_{t}\boldsymbol{ \phi }$ term in (\oldref{eqS2regs}). For the \emph{`correctly specified}\textit{'} Stage 2 model in (\oldref{S2full0}), $\boldsymbol{\mathcal{X}}_{t}$ in (\oldref{XX}) is replaced by the $(1\times 3)$ vector $[\hat{\tilde{y}}_{t-1|T},~\hat{ \tilde{y}}_{t-2|T},~(r_{t-1}+r_{t-2}-4\{\hat{g}_{t-1|T}+\hat{g} _{t-2|T}\})/2] $.

In the second sequence, labelled `constant $\boldsymbol{\phi }$ ' in (ref) and drawn as a blue line, the observable counterpart to $GY_{t}$ is computed only once outside the structural break regression loop, with the dummy variable regression performed without the extra $\boldsymbol{\mathcal{X}}_{t}\boldsymbol{\phi }$ term in (\oldref{eqS2regs}) , ie., it is computed in its `original' form as given in (\oldref{Zt}). \footnote{ Note that \cites{stock.watson:1998} MUE look-up table values for $\lambda $ were constructed by simulation with the structural break test testing the unconditional mean of the $GY_{t}$ series for a break, without any other variables being included in the regression. This form of the structural break regression is thus compatible with \cites{stock.watson:1998} look-up table values.} More specifically, for the `misspecified ' and `correctly specified\textit{' }Stage 2 models, the observable counterparts to the $GY_{t}$ series are constructed as:

respectively. The $\hat{a}_{y,1},\hat{a}_{y,2},\hat{a}_{r},\hat{a}_{g},$ and $\hat{a}_{0}$ coefficients are the (full sample) estimates reported in columns 2 and 4 of (ref) under the headings `Replicated' and `MLE$(\sigma _{g}).\mathcal{M}_{0}$', with the corresponding latent state estimates from the respective models.\footnote{ For instance, $\hat{g}_{t-1|T}$ in (\oldref{GY_HLW}) is the Kalman Smoothed estimate of trend growth from \cites{holston.etal:2017} `misspecified ' Stage 2 model, while trend growth $\hat{g}_{t-1|T}$ in (\oldref{GYcorr1}) is the corresponding estimate from the `correctly specified ' Stage 2 model.}

As can be seen from (ref), the $\{F(\tau )\}_{\tau =\tau _{0}}^{\tau _{1}}$ sequences from the `correctly specified' Stage 2 models shown in the bottom panel are not only smaller overall, but they are nearly unaffected by \cites{holston.etal:2017} approach to `re-estimate' the parameters in the structural break loop. Both, the `constant $\boldsymbol{\phi }$' and the `time varying $\boldsymbol{\phi }$\textit{'} versions generate $ \{F(\tau )\}_{\tau =\tau _{0}}^{\tau _{1}}$ sequences that are overall very similar, with their maximum values being around 4.5. For the \emph{ `misspecified}\textit{'} Stage 2 model shown in the top panel, this is not the case. The variation as well as the magnitude of $\{F(\tau )\}_{\tau =\tau _{0}}^{\tau _{1}}$ from the \emph{`time varying} $\boldsymbol{\phi }$ \textit{'} and \emph{`constant} $\boldsymbol{\phi }$\textit{' } implementations are vastly different, with the former having a much higher mean and maximum value.

These large differences in the $\{F(\tau )\}_{\tau =\tau _{0}}^{\tau _{1}}$ sequences from the `misspecified' Stage 2 models also lead to very different estimates of $\lambda _{z}$. This can be seen from (ref), which shows the resulting $\lambda _{z}$ estimates in the top part with the corresponding $L$, MW, EW, and QLR structural break test statistics in the bottom part. (ref) is arranged further into a left and a right column block, referring to the `time varying $\boldsymbol{\phi }$\textit{' } and the \emph{`constant} $\boldsymbol{\phi }$\textit{' }MUE implementations for the three different models reported in (\oldref{tab:Stage2}). `Replicated' refers to the baseline replicated results, `MLE$(\sigma _{g})$' corresponds to the \emph{`misspecified}\textit{'} Stage 2 model but with $\sigma _{g}$ estimated by MLE, and `MLE$(\sigma _{g}).\mathcal{M}_{0}$' is from the \emph{ `correctly specified}\textit{'} Stage 2 model with $\sigma _{g}$ again estimated by MLE. The `HLW.R-File' column lists the results from \cites{holston.etal:2017} R-Code. Note that holston.etal:2017 do not report estimates based on \cites{nyblom:1989} $L$ statistic. The entries in the $L$ rows in (ref) under `HLW.R-File' thus simply list `---'. 90% confidence intervals for $\lambda _{z}$ and $p-$values for the structural break tests are reported in square and round brackets, respectively. \footnote{ As in the replication of \cites{stock.watson:1998} results reported in (\oldref{tab:sw98_T4}), these were again obtained from their GAUSS files.}

Consistent with the visual findings from (ref), the structural break statistics from the `misspecified' Stage 2 model shown under the 'Replicated' heading for the `time varying $\boldsymbol{ \phi }$' and `constant $\boldsymbol{\phi }$' settings are very different. The \textrm{MW}, \textrm{EW}, and \textrm{QLR} statistics are approximately 4 to 5 times larger under the \emph{`time varying} $\boldsymbol{\phi }$\textit{' }setting than under the \emph{ `constant} $\boldsymbol{\phi }$\textit{'} scenario. Because \cites{nyblom:1989} $L$ statistic is constructed as the scaled cumulative sum of the demeaned `$GY_{t}$' series and thus does not require the partitioning of data, creation of dummy variables, or looping through potential break dates, it is not affected by this choice, yielding the same test statistic of about $0.05$ under both settings.

Under the `time varying $\boldsymbol{\phi }$' setting, the MW, EW, and QLR statistics and \cites{nyblom:1989} $L$ statistic generate vastly different $\lambda _{z}$ estimates. \cites{nyblom:1989} $L$ statistic is highly insignificant with a $p-$value of $0.87$, resulting in a $\lambda _{z}$ estimate of exactly $0$ ( \cites{nyblom:1989} $L$ statistic is less than $0.118$, the smallest value in \cites{stock.watson:1998} look-up Table 3 which corresponds to $\lambda =0$). The MW, \textrm{ EW}, and \textrm{QLR} structural break statistics on the other hand are either weakly significant or marginally insignificant, with $p-$values between $0.045$ and $0.13$. These borderline significant structural break statistics generate sizable $\lambda _{z}$ point estimates between $0.025$ and $0.034$. The resulting $90\%$ confidence intervals for $\lambda _{z}$ are, nonetheless, rather wide with $0$ as the lower bound, suggesting that these point estimates are not significantly different from zero.\footnote{ Given the earlier discussion in \Sref{subsec:MUE} and the ARE results in Table 2 of stock.watson:1998, we know that MUE\ can be a very inefficient estimator.} Under the \emph{`constant} $\boldsymbol{\phi }$ \textit{'} setting, the four structural break statistics and the resulting $ \lambda _{z}$ estimates tell a consistent story (see the `Replicated' heading in the right column block). All structural break statistics are highly insignificant, with their respective $\lambda _{z}$ point estimates being equal to zero.

For the `correctly specified' Stage 2 models shown under the headings `MLE$(\sigma _{g}).\mathcal{M}_{0}$' in (ref), the `time varying $\boldsymbol{\phi }$' and the `constant $\boldsymbol{\phi }$' estimates of $ \lambda _{z}$ reflect the visual similarity of the $\{F(\tau )\}_{\tau =\tau _{0}}^{\tau _{1}}$ sequences shown in the bottom panel of (ref) . The $\lambda _{z}$ point estimates are of the same order of magnitude, very close to zero (they are exactly equal to zero for \cites{nyblom:1989} $ L $ statistic and \textrm{MW}\ under the \emph{`constant} $\boldsymbol{\phi } $\textit{'} setting), and most importantly, substantially smaller than those constructed from \cites{holston.etal:2017} \textit{`}\emph{misspecified}\textit{'} Stage 2 model.\footnote{ In (ref) in the \hyperref[appendix]{Appendix}, I present these Stage 2 MUE results for data that was updated to 2019:Q2. The conclusion is the same.}

What is causing this large difference in the $\{F(\tau )\}_{\tau =\tau _{0}}^{\tau _{1}}$ sequences between the\ `misspecified' and `correctly specified' Stage 2 models in the `time varying $\boldsymbol{\phi }$' setting? There are two components. First, the Kalman Smoothed estimates of the output gap (cycle) $\hat{\tilde{y }}_{t|T}\ $and of (annualized) trend growth $\hat{g}_{t|T}$ can be quite different from these two models, despite the parameter estimates and values of the log-likelihoods being very similar. This difference is more pronounced for the cycle estimate $\hat{\tilde{y}}_{t|T}$, particulary towards the end of the sample period than for the trend growth estimate $ \hat{g}_{t|T}$ (see (ref) in the \hyperref[appendix]{ Appendix} which shows a comparison of $\hat{\tilde{y}}_{t|T}\ $and $\hat{g} _{t|T}$ from the\ \emph{`misspecified}\textit{'} and \emph{`correctly specified}\textit{'} Stage 2 models).

Second, the parameter restriction $(a_{g}+4a_{r})$ on the relationship between the real rate and trend growth matters. More specifically, when conditioning on $\boldsymbol{\mathcal{X}}_{t}$ in (\oldref{eqS2regs}), it is the restriction $(r_{t-1}-4\hat{g}_{t-1|T})$ in $\boldsymbol{\mathcal{X}}_{t}$ that makes the largest difference to the $\{F(\tau )\}_{\tau =\tau _{0}}^{\tau _{1}}$ sequence. To see this, I\ show plots of the $\{F(\tau )\}_{\tau =\tau _{0}}^{\tau _{1}}$ sequences from various $\boldsymbol{ \mathcal{X}}_{t}$ constructs corresponding to the different Stage 2 model specifications in (ref) in the \hyperref[appendix]{Appendix} . I use the `correctly specified' Stage 2 model's $\{\hat{ \tilde{y}}_{t-i|T}\}_{i=1}^{2}$ and $\hat{g}_{t-1|T}$ estimates to form three sets of $\boldsymbol{\mathcal{X}}_{t}$ vectors for the dummy variable regressions in (\oldref{eqS2regs}). These are:\bsq

\esq and are labelled accordingly in (ref) (the preceding `MLE$(\sigma _{g}).\mathcal{M}_{0}$' signifies that these were constructed using the $\{\hat{\tilde{y}}_{t-i|T}\}_{i=1}^{2}$ and $\hat{g}_{t-1|T}$ estimates from the `correctly specified' Stage 2 model). The corresponding $\mathcal{Y}_{t}$ dependent variable for these structural break regressions also uses the `correctly specified' Stage 2 model's output gap estimate $\hat{\tilde{y}}_{t|T}$. The $\{F(\tau )\}_{\tau =\tau _{0}}^{\tau _{1}}$ sequences from \cites{holston.etal:2017} `misspecified' and the `correctly specified\textit{'} Stage 2 models are superimposed as reference values and are denoted by `HLW' and `MLE$(\sigma _{g}).\mathcal{M}_{0}$'.

The plot corresponding to (\oldref{c1}) (orange dashed line in (ref)) shows a rather small difference relative to the `HLW' benchmark (blue solid line). Thus, exchanging $\{\hat{\tilde{y}} _{t-i|T}\}_{i=1}^{2}$ and $\hat{g}_{t-1|T}$ from \cites{holston.etal:2017} `misspecified' Stage 2 model for those from the `correctly specified' one only has a small impact on the $\{F(\tau )\}_{\tau =\tau _{0}}^{\tau _{1}}$ sequence and is most visible over the 1994 to 2000 period. Dropping the second lag in $r_{t}$ from $\boldsymbol{ \mathcal{X}}_{t}$ in (\oldref{c2}) (see the cyan dotted line in (ref)) also has only a small impact on the $\{F(\tau )\}_{\tau =\tau _{0}}^{\tau _{1}}$ sequence. The biggest effect on $\{F(\tau )\}_{\tau =\tau _{0}}^{\tau _{1}}$ has the restriction $(r_{t-1}-4\hat{g}_{t-1|T})$ as imposed in (\oldref{c3}) (green dashed-dotted line (ref)). This is evident from the near overlapping with the red solid line corresponding to the correctly specified' Stage 2 model's $\{F(\tau )\}_{\tau =\tau _{0}}^{\tau _{1}}$ sequence. Recall that the only difference between these two is that an extra lag of $(r_{t-1}-4\hat{g}_{t-1|T})$ is added to $\boldsymbol{\mathcal{X}}_{t}$, and that these enter as an average, viz, $\boldsymbol{\mathcal{X}}_{t}=[\hat{\tilde{y}}_{t-1|T},~\hat{\tilde{y}} _{t-2|T},~(r_{t-1}+r_{t-2}-4\{\hat{g}_{t-1|T}+\hat{g}_{t-2|T}\})/2]$.

What does \cites{holston.etal:2017} Stage 2 MUE procedure recover?

\cites{holston.etal:2017} Stage 2 MUE procedure implemented on the `misspecified' Stage 2 model leads to spuriously large estimates of $\lambda _{z}$ when the true value is zero. To show this, I\ perform two simple simulation experiments.

In the first experiment, I\ simulate data from the full structural model in (\oldref{eq:hlw}) using the Stage 3 parameter estimates of holston.etal:2017 reported in column one of (ref) as the true values that generate the data, but with `other factor' $z_{t}$ set to zero for all $t$. The natural rate $r_{t}^{\ast }$ in the output gap equation in (\oldref{IS}) is thus solely determined by (annualized) trend growth, that is, $r_{t}^{\ast }=4g_{t}$, which implies that $\lambda _{z}$ is zero in the simulated data.\footnote{ To implement the simulations from the full Stage 3 model, I\ need to define a process for the exogenously determined interest rate in \cites{holston.etal:2017} model. For simplicity, I estimate a parsimonious, but well fitting, ARMA($2,1$) model for the real interest rate series, and then use the ARMA($2,1$) coefficients to generate a sequence of 229 simulated observations for $r_{t}$. Recall that holston.etal:2017 use data from 1960:Q1, where the first 4 quarters are used for initialisation of the state vector, so that in total $4+225=T$ observations are available. The remaining series are simulated from the Stage 3 model given in (\oldref{eq:hlw}). To get a realistic simulation path from the Stage 3 model, I\ initialize the first four data points for the simulated inflation series at their observed empirical values. For the $y_{t}^{\ast }$ series, the HP-filter based trend estimates of GDP (also utilized in the initialisation of the State vector in Stage 1) are used to set the first four observations. The cycle variable $ \tilde{y}_{t}$ is initialized at zero, while trend growth $g_{t}$ is initialized at $0.75$, which corresponds to an annualized rate of $3$ percent. In the analysis that requires a simulated path of `other factor' $z_{t}$, ie., when the natural rate is generated from $r_{t}^{\ast }=4g_{t}+z_{t}$, the first four entries in $z_{t}$ are initialized at zero. A total of $S=1000$ sequences are simulated with a total sample size of $229$ observations, where the first four entries are discarded in later analysis.} I then implement \cites{holston.etal:2017} Stage 2 MUE procedure on the simulated data following steps (I\hsp[.3]) to (IV\hsp[.3]) outlined in \Sref{sec:MUE2}\ above to yield a sequence of $S=1000$ estimates of $\lambda _{z}$ $\left( \{\hat{\lambda}_{z}^{s}\}_{s=1}^{S}\right) $.

I use two different scenarios for $\boldsymbol{\theta }_{2}$ in the Kalman Smoother recursions described in Step $(I)$ to extract the latent cycle as well as trend growth series needed for the construction of $\mathcal{Y}_{t}$ and $\boldsymbol{\mathcal{X}}_{t}$ in the dummy variable regression in (\oldref{eqS2regs}). The first scenario simply takes \cites{holston.etal:2017} empirical Stage 2 estimate $\skew{0}\boldsymbol{\hat{\theta}}_{2}$ as reported in column one of (ref), and keeps these values fixed for all 1000 generated data sequences when applying the Kalman Smoother. In the second scenario, I re-estimate the Stage 2 parameters for each simulated sequence to obtain new estimates $\skew{0}\boldsymbol{\hat{\theta}} _{2}^{s},\forall s=1,\ldots ,S$. I then apply the Kalman Smoother using these estimates to generate the $\mathcal{Y}_{t}$ and $\boldsymbol{\mathcal{X }}_{t}$ sequences for the regression in (\oldref{eqS2regs}).

Finally, I repeat the above computations on data that were generated from the full model in (\oldref{eq:hlw}) with the natural rate of interest determined by both factors, namely, $r_{t}^{\ast }=4g_{t}+z_{t}$, where $z_{t}$ was simulated as a pure random walk. The standard deviation of $z_{t}$ was set at the implied value from the Stage 2 estimate of $\lambda _{z}$ and the Stage 3 estimates of $\sigma _{\tilde{y}}$ and $a_{r}$, ie., at $\sigma _{z}=\lambda _{z}\sigma _{\tilde{y}}/a_{r}\approx 0.15$ (see row $\sigma _{z} $ (implied) of column one in (ref)). The objective here is to provide a comparison of the magnitudes of the $\lambda _{z}$ estimates that are obtained when implementing \cites{holston.etal:2017} Stage 2 MUE procedure on data that were generate with and without `other factor' $z_{t}$ in the natural rate.

In (ref), summary statistics of $\hat{\lambda} _{z}^{s}$ from the two different data generating processes (DGPs) are reported. The left column block shows results for the two different DGPs when the Stage 2 parameter vector $\boldsymbol{\theta }_{2}$ is held fixed at the estimates reported in column one of (ref). The right column block shows corresponding results when $\boldsymbol{\theta }_{2}$ is re-estimated for each simulated data series. The summary statistics are the minimum, maximum, standard deviation, mean, and median of $\hat{\lambda} _{z}^{s}$, as well as the relative frequency of obtaining a value larger than the empirical point estimate of holston.etal:2017.\ This point estimate and the corresponding relative frequency are denoted by $\hat{ \lambda}_{z}^{\mathrm{HLW}}$ and $\Pr (\hat{\lambda}_{z}^{s}>\hat{\lambda} _{z}^{\mathrm{HLW}})$, respectively. To complement the summary statistics in (ref), histograms of $\hat{\lambda}_{z}^{s}$ are shown in (ref) to provide visual information about its sampling distribution.

From the summary statistics in (ref) as well as the histograms in (ref) we can see how similar the $\hat{ \lambda}_{z}^{s}$ coefficients from these two different DGPs are. For instance, when the data were simulated without `other factor' $ z_{t} $ (ie., $\lambda _{z}=0$), the sample mean of $\hat{\lambda}_{z}^{s}$ is $0.028842$. When the data were generated from the full model with $ r_{t}^{\ast }=4g_{t}+z_{t}$, the sample mean of $\hat{\lambda}_{z}^{s}$ is only $6.53\%$ higher at $0.030726$. Similarly, the relative frequencies $\Pr (\hat{\lambda}_{z}^{s}>\hat{\lambda}_{z}^{\mathrm{HLW}})$ for these two DGPs are $45.70\%$ and $49\%$, respectively. The inclusion of `other factor' $z_{t}$ in the DGP\ of the natural rate thus results in only a $3.3$ percentage points higher $\Pr (\hat{\lambda}_{z}^{s}>\hat{\lambda}_{z}^{ \mathrm{HLW}})$.\footnote{ When the Stage 2 parameter vector $\boldsymbol{\theta }_{2}$ is re-estimated for each simulated sequence shown in the right column block in (ref), the sample means as well as the relative frequency $ \Pr (\hat{\lambda}_{z}^{s}>\hat{\lambda}_{z}^{\mathrm{HLW}})$ are somewhat lower at $0.025103$ and $0.027462$, and 33.90% and 39.30%, respectively.} The histograms in (ref) paint the same overall picture. As can be seen, the Stage 2 MUE implementation has difficulties to discriminate between these two DGPs. Moreover, it seems that it is \cites{holston.etal:2017} procedure itself that leads to the spuriously amplified estimates of $ \lambda _{z}$, regardless of the data.

In a second experiment I\ simulate DGPs from entirely unrelated univariate ARMA\ processes of the individual components of the $\mathcal{Y}_{t}$ and $ \boldsymbol{\mathcal{X}}_{t}$ series needed for the regressions in (\oldref{eqS2regs}). To match the time series properties of the $\mathcal{Y}_{t}$ and $\boldsymbol{\mathcal{X}}_{t}$ elements given in (\oldref{YY}) and (\oldref{XX}), I fit simple low-order ARMA\ models to $\hat{\tilde{y}}_{t|T},$ $r_{t}$ and $ \hat{g}_{t|T}$, and then use these ARMA\ estimates to simulate artificial data.\footnote{ I use 4 different time series processes for $\hat{g}_{t|T}$ in these simulations. Complete details of the simulation design are given in \hyperref[sec:AS4]{Section A.4} of the \hyperref[appendix]{Appendix}.} Finally, I apply \cites{holston.etal:2017} Stage 2 MUE procedure to the simulated data as before, nevertheless starting from Step (II\hsp[.3] ), and thereby skipping the Kalman Smoother step. The full results from the second experiment are reported in (ref) and (ref) in the \hyperref[appendix]{Appendix}. These yield magnitudes of $\hat{\lambda}_{z}^{s}$ that are similar to those from the first simulation experiment, with mean estimates being between $0.026117$ and $0.031798$, and relative frequencies corresponding to $\Pr (\hat{\lambda} _{z}^{s}>\hat{\lambda}_{z}^{\mathrm{HLW}})$ being between $38.40\%$ and $ 49.80\%$.

Stage 3 Model

The analysis so far has demonstrated that the ratios of interest $\lambda _{g}=\sigma _{g}/\sigma _{_{y^{\ast }}}$ and $\lambda _{z}=a_{r}\sigma _{z}/\sigma _{\tilde{y}}$ required for the estimation of the full structural model in (\oldref{eq:hlw}) cannot be recovered from \cites{holston.etal:2017} MUE\ procedure implemented in Stages 1 and 2. Moreover, since their procedure is based on the `misspecified' Stage 2 model in (\oldref{eq:stag2}), it results in a substantially larger estimate of $\lambda _{z}$ than when implemented on the `correctly specified' Stage 2 model in (\oldref{S2full0}). This substantially larger estimate of $ \lambda _{z}$ in turn leads to a greatly amplified and strongly downward trending `other factor' $z_{t}$. To show the impact of this on \cites{holston.etal:2017} estimate of the natural rate of interest, I initially report parameter estimates of the full Stage 3 model in (ref), followed by plots of filtered estimates of the natural rate $ r_{t}^{\ast }$, trend growth $g_{t}$, `other factor' $z_{t}$, and the output gap (cycle) variable $\tilde{y}_{t}$ in (ref). \footnote{ Smoothed estimates are shown in (ref). In \hyperref[sec:AS3]{ Section A.3} in the \hyperref[appendix]{Appendix}, the expansion of the system matrices are reported as for the earlier Stage 1 and Stage 2 models. These are in line with the full model reported in (\oldref{eq:hlw}). As before, the state vector $\boldsymbol{\xi }_{t}$ is initialized using the same procedure as outlined in (\oldref{eq:P00S1a}) and \fnref{fn:1}, with the numerical values of $\boldsymbol{\xi }_{00}$ and $\mathbf{P}_{00}$ given in (\oldref{AS3:xi00}) and (\oldref{AS3:P00}).}

Given estimates of the ratios $\lambda _{g}=\sigma _{g}/\sigma _{_{y^{\ast }}}$ and $\lambda _{z}=a_{r}\sigma _{z}/\sigma _{\tilde{y}}$ from the previous two stages, the vector of Stage 3 parameters to be computed by MLE is:

In (ref), estimates of $\boldsymbol{\theta }_{3}$ are presented following the same format as in (ref) and (ref) previously. Since I also estimate $\sigma _{g}$ and $\sigma _{z} $ directly together with the other parameters by MLE without using the Stage 1 and Stage 2 estimates of $\lambda _{g}$ and $\lambda _{z}$, additional rows are inserted, with the values in brackets denoting implied estimates. The first two columns in (ref) show estimates of $ \boldsymbol{\theta }_{3}$ obtained from running \cites{holston.etal:2017} R-Code and my replication. The third and fourth columns (under headings `MLE( $\sigma _{g}|\hat{\lambda}_{z}^{\mathrm{HLW}}$)' and `MLE($\sigma _{g}|\lambda _{z}^{\mathcal{M}_{0}}$)', respectively) report estimates when $ \sigma _{g}$ is estimated freely by MLE, while $\lambda _{z}$ is held fixed at either $\hat{\lambda}_{z}^{\mathrm{HLW}}=0.030217$ obtained from \cites{holston.etal:2017} `misspecified' Stage 2 model under their `time varying $\boldsymbol{\phi }$' approach, or at $ \hat{\lambda}_{z}^{\mathcal{M}_{0}}=0.000754$ computed from the `correctly specified' Stage 2 model in (\oldref{S2full0}) with \emph{ `constant} $\boldsymbol{\phi }$\textit{'}. The last column of (ref) under heading `MLE($\sigma _{g},\sigma _{z}$)' lists the estimates of $\boldsymbol{\theta }_{3}$ when $\sigma _{g}$ and $\sigma _{z}$ are computed directly by MLE, with the implied values of $\lambda _{g}$ and $ \lambda _{z}$ reported in brackets.

The Stage 3 results in (ref) can be summarized as follows.\ The MLE of $\sigma _{g}$ does not `pile-up' at zero and is again approximately $50\%$ larger than the estimate implied by the Stage 1 MUE\ of $\lambda _{g}$. That is, $\hat{\sigma}_{g}\approx $ $0.045$ in the last three columns of (ref), and thus very similar in size to the Stage 2 estimates of $0.044$ and $0.045$ shown in the last two columns of (ref). Computing $\sigma _{z}$ directly by MLE leads to a point estimate that shrinks numerically to zero, while the estimates of the other parameters remain largely unchanged. Notice again that the log-likelihood values of the last three models in (ref) are very similar, ie., between $-514.8307$ and $-514.2899$. Yet, the corresponding estimates of $\sigma _{z}$ are either very small at $0$ or comparatively large at $0.1371$ when implied from the `misspecified ' Stage 2 model's $\hat{\lambda}_{z}^{\mathrm{HLW}}$ estimate. The $ \hat{\sigma}_{z}$ coefficient from the `correctly specified' Stage 2 model is $0.0037$ and thereby nearly 40 times smaller than from the `misspecified\textit{'} Stage 2 model.

The findings from (ref) are mirrored in the filtered estimates of $r_{t}^{\ast }$, $g_{t}$, $z_{t}$ and $\tilde{y}_{t}$ plotted in (ref). The `MLE($\sigma _{g}|\lambda _{z}^{\mathcal{M} _{0}} $)' and `MLE($\sigma _{g},\sigma _{z}$)' estimates are visually indistinguishable. Unsurprisingly, out of the four estimates, `other factor' $z_{t}$ is overall most strongly affected by the two different $ \lambda _{z}$ values that are conditioned upon, showing either\ vary large variability and a pronounced downward trend in $z_{t}$, or being close to zero with very little variation (see panel (c) in (ref)). The effect on the estimate of the natural rate is largest in the immediate aftermath of the global financial crisis, namely, from 2010 onwards. Interestingly, the output gap estimates shown in panel (d) of (ref) are quite similar, with the largest divergence occurring after 2012. The three trend growth estimates in panel (b) of (ref) which estimate $\sigma _{g}$ directly by MLE are visually indistinguishable, despite having very different ${\sigma}_{z}$ values, namely, between $0$ and $0.1371$ (see the lines corresponding to `MLE($\sigma _{g}|\lambda _{z}^{ \mathrm{HLW}}$)', `MLE($\sigma _{g}|\lambda _{z}^{\mathcal{M}_{0}}$)' and `MLE($\sigma _{g},\sigma _{z}$)'). Trend growth estimated from \cites{holston.etal:2017} \ Stage 1 MUE of $\lambda _{g}$ is noticeably larger from 2009 to 2014. In comparison to the plots shown in panel (c) of (ref), the drop in all four trend growth estimates following the financial crisis seems exaggerated. The pure backward looking nature of the Kalman Filtered $g_{t}$ series exacerbates the effect of the decline in GDP during the financial crisis on trend growth estimates after the crisis.

A final point I would like to make here --- and without the intention to engage in repetitive and unnecessary discussion --- is that extending the sample period to 2019:Q2 produces the interesting empirical result that estimating $\sigma _{z}$ directly by MLE does not lead to any `pile-up' at zero problems. Moreover, the ML estimate of $\sigma _{z}$ is very similar to the one implied from the `correctly specified' Stage 2 model's $\lambda _{z}$, and thereby again in stark contrast to the oversized estimate obtained from \cites{holston.etal:2017} `misspecified' Stage 2 model's $\lambda _{z}$.\footnote{ These estimation results using data up to 2019:Q2 together with corresponding plots of filtered (and smoothed) estimates are reported in (ref), (ref) and (ref) in \hyperref[sec:AS3]{Section A.3} of the \hyperref[appendix]{Appendix}.} Even so, despite the fact that the point estimate of $\sigma _{z}$ does not shrink to zero, it is highly insignificant, which suggests that there is little evidence in the data of `other factor' $z_{t}$ being relevant for the model.

\stepcounter{section} \Osubsection*{A.\arabic{section}. {Other issues}} \addcontentsline{toc}{subsection}{\currentname}

There are other issues with \cites{holston.etal:2017} structural model in (\oldref{eq:hlw}) that make it unsuitable for policy analysis. For instance, the interest rate $i_{t}$ is included as an exogenous variable, so that the model essentially tries to find the best fitting natural rate $r_{t}^{\ast }$ for it. With $r_{t}^{\ast }=4g_{t}+z_{t} $, and `other factor' $z_{t}$ the `free' variable due to $g_{t}$ being driven by GDP, $z_{t}$ effectively matches the `leftover' movements in the interest rate to make it compatible with trend growth in the model. Since the central bank has full control over the (fed funds) interest rate, it can set $i_{t}$ to any desired level and the model will produce a natural rate through `other factor' $z_{t}$ that will match it. Also, there is nothing in the structural model of (\oldref{eq:hlw}) that makes the system stable. For the output gap relation in (\oldref{IS}) to be stationary, the real rate cycle $r_{t}-r_{t}^{\ast }=(i_{t}-\pi _{t}^{e})-(4g_{t}+z_{t})$ must be $I(0)$, yet there is no co-integrating relation imposed anywhere in the system to ensure that this holds in the model.\footnote{ This insight is not new and has been discussed in, for instance, pagan.wickens:2019 (see pages $21-23$).} When trying to simulate from such a model, with $\pi _{t}$ being integrated of order 1, the simulated paths of the real rate $r_{t}=i_{t}-\pi _{t}^{e}$ can frequently diverge to very large values, even with samples of size $T=229$ observations, which is the empirical sample size.

A broader concern for policy analysis is the fact that the filtered estimates of the state vector $\boldsymbol{\xi }_{t}$ will be (weighted combinations of the) one-sided moving averages of the three observed variables that enter the state-space model; namely, $i_{t}$, $y_{t}$, and $ \pi _{t}$.\footnote{ Smoothed estimates will be (weighted combinations of the) two-sided moving averages of the observables. See also durbin.koopman:2012, who write to this on page 104: "It follows that these conditional means are weighted sums of past (filtering), of past and present (contemporaneous filtering) and of all (smoothing) observations. It is of interest to study these weights to gain a better understanding of the properties of the estimators as is argued in Koopman and Harvey (2003). ... . In effect, the weights can be regarded as what are known as kernel functions in the field of nonparametric regression; ... ."} This can be seen by writing out the Kalman Filtered estimate of the state vector as:\footnote{ I again follow the notation in hamilton:1994, see pages 394-395, with the matrices $\mathbf{A}$ and $\mathbf{H}$ however not transposed to be consistent with my earlier notation.}

where $\boldsymbol{\Psi }_{i}=\prod_{n=0}^{i-1}\mathbf{\Phi }_{t-n},\forall i=1,2,\ldots ,$ $\boldsymbol{\Psi }_{0}=\mathbf{I,}$ $\mathbf{I}$ is the identity matrix, $\skew{0}\boldsymbol{\hat{\xi}}_{t|t-1}=\mathbf{F}\skew{0} \boldsymbol{\hat{\xi}}_{t-1|t-1}$ is the predicted state vector, $ \boldsymbol{\xi }_{0|0}$ is the prior mean, $\mathbf{P}_{t|t-1}=\mathbf{FP} _{t-1|t-1}\mathbf{F}+\mathbf{Q}$ is the predicted state variance, $ \boldsymbol{\omega }_{ti}=\boldsymbol{\Psi }_{i}\mathbf{G}_{t-i}$ is a time varying weight matrix, and $\mathbf{\bar{y}}_{t}$ consists of the observed variables $y_{t}$, $\pi _{t}$, and $i_{t}$.\footnote{ To understand what is driving the downward trend in `other factor' $ z_{t}$ since the early 2000s in the model, one could examine the weight matrix $\boldsymbol{\omega }_{ti}$ in (\oldref{xirec}) more closely to see how it interacts with the observable vector $\mathbf{\bar{y}}_{t}=\mathbf{y}_{t}- \mathbf{Ax}_{t}=[a(L)y_{t}-a_{r}(L)r_{t};~b_{\pi }(L)\pi _{t}-b_{y}y_{t}]^{\prime }$, where $b_{\pi }(L)=1-b_{\pi }L-\frac{1}{3} (1-b_{\pi })(L^{2}+L^{3}+L^{4})$ is the lag polynomial capturing the dynamics of inflation. Alternatively, the steady-state $\mathbf{P}$ matrix could be computed recursively as in equation 13.5.3 in hamilton:1994 to replace $\mathbf{P}_{t|t-1}$ in the recursions for $\skew{0}\boldsymbol{ \hat{\xi}}_{t|t}$. The relation in (\oldref{xirec}) would then yield $\skew{0} \boldsymbol{\hat{\xi}}_{t|t}=\Phi ^{t}\boldsymbol{\xi }_{0|0}+ \sum_{i=0}^{t-1}\mathbf{\Phi }^{i}\mathbf{G\bar{y}}_{t-i}$, where $\mathbf{ \Phi =(I-GH)F}$ and $\mathbf{G=PH}^{\prime }(\mathbf{HP}^{\prime }\mathbf{H} ^{\prime }+\mathbf{R})^{-1}$ would be the steady-state analogue to $\mathbf{ \Phi }_{t}$ and $\mathbf{G}_{t}$, with $\mathbf{P}_{t|t-1}$ replaced by $ \mathbf{P}$ from the steady-state $\mathbf{P}$ matrix.}

This creates the following two issues. First, since the nominal interest rate $i_{t}$ is directly controlled by the central bank, and the natural rate is constructed from the filtered estimate of state vector $\boldsymbol{ \xi }_{t}$, which itself is computed as a moving average of $i_{t}$ (and the other observable variables), a circular relationship can be seen to evolve. Any central bank induced change in the policy rate $i_{t}$ is mechanically transferred to the natural rate $r_{t}^{\ast }$ via the Kalman Filtered estimate of the state vector $\skew{0}\boldsymbol{\hat{\xi}}_{t|t}$ in (\oldref{xirec}). A confounding effect between $r_{t}^{\ast }$ and $i_{t}$ will arise, making it impossible to answer questions of interest such as: \textquotedblleft Is the natural rate low because $i_{t}$\ is low, or is $i_{t}$\ low because the natural rate is low?\textquotedblright with this model, as one will follow as a direct consequence from the other.

Second, because of the one-sided moving average nature of the Kalman Filtered estimates of the state vector, any outliers, structural breaks or otherwise `extreme' observations at the beginning (or end) of the sample period can have a strong impact on these filtered estimates. For the (two-sided) hodrick.prescott:1997 filter, such problems (and other ones) are well known and have been discussed extensively in the literature before.\footnote{ There exists a large literature on the HP filter and its problems (one of the more recent papers is by hamilton:2018), and it is not the goal to review or list them here. However, the study by phillips.jin:2015 is interesting to single out, in particular the introduction section on pages 2 to 9, as it highlights the recent public debates by James Bullard, Paul Krugman, Tim Duy and others on the use (and misuse) of the HP filter for the construction of output gaps for policy analysis. phillips.jin:2015 show also that the HP filter fails to recover the underlying trend asymptotically in models with breaks (see section 4 in their paper), and they further propose alternative filtering/smoothing methods. In an earlier study, schlicht:2008 describes how to deal with structural breaks and missing data.} However, (one-sided) Kalman Filter based estimates will also be affected. This can be easily demonstrated here by re-estimating the model using four different starting dates, while keeping the end of the sample period the same at 2019:Q2. In (ref) I show filtered estimates of $r_{t}^{\ast }$, $g_{t}$, $z_{t}$ and $\tilde{y}_{t}$ for the starting dates 1967:Q1, 1972:Q1, 1952:Q2 and 1947:Q1 (smoothed estimates are shown in (ref)), together with \cites{holston.etal:2017} estimates using 1961:Q1 as the starting date. \footnote{ In all computations, I use \cites{holston.etal:2017} R-Code and follow exactly their three stage procedure as before to estimate the factors of interest.}

Why are these starting dates chosen? The period following the April 1960 to February 1961 recession was marked by temporarily unusually (and perhaps misleadingly) high GDP\ growth, yielding an annualized mean of $6.07\%$ (median $6.47\%$), with a low standard deviation of $2.67\%$ from 1961:Q2 to 1966:Q1 (see panel (b) of (ref)). Having such `excessive' growth at the beginning of the sample period has an unduly strong impact on the filtered (less so on the smoothed) estimate of trend growth $g_{t}$ in the model. Since both $g_{t}$ and $z_{t}$ enter the natural rate, this affects the estimate of $r_{t}^{\ast }$. To illustrate the sensitivity of these estimates to this time period, I\ re-estimate the model with data starting 6 years later in 1967:Q1. Also, \cites{holston.etal:2017} Euro Area estimates of $r_{t}^{\ast }$ are negative from around 2013 onwards (see the bottom panel of Figure 3 on page S63 of their paper).\footnote{ This negative estimate in $r_{t}^{\ast }$ is driven by an excessively large and volatile estimate of `other factor' $z_{t}$. Some commentators have attributed the larger decline in the natural rate to a stronger manifestation of `secular stagnation' in the Euro Area than in the U.S.} To show that we can get the same negative estimates of $r_{t}^{\ast }$ for the U.S., I re-estimate the model with data starting in 1972:Q1 to match the sample period of the Euro Area in holston.etal:2017. Lastly, I\ extend \cites{holston.etal:2017} data back to 1947:Q1 to have estimates from a very long sample, using total PCE\ inflation prior to 1959:Q2 in place of Core PCE\ inflation and the Federal Reserve Bank of New York discount rate from 1965:Q1 back to 1947:Q1 as a proxy for the Federal Funds rate as was done in laubach.williams:2003. \footnote{ Note that from the quarterly CORE PCE data it will be possible to construct annualized inflation only from 1947:Q2 onwards. To have an inflation data point for 1947:Q1, annual core PCE\ data (BEA\ Series ID: DPCCRG3A086NBEA) that extends back to 1929 was interpolated to a quarterly frequency and subsequently used to compute (annualized) quarterly inflation data for 1947:Q1. Since \cites{holston.etal:2017} R-Code requires 4 quarters of GDP\ data prior to 1947:Q1 as initial values, annual GDP\ (BEA Series ID: GDPCA) was interpolated to quarterly data for the period 1946:Q1 1946:Q4.} Since inflation was rather volatile from 1947 to 1952, I also re-estimate the model with data beginning in 1952:Q2 to exclude this volatile inflation period from the sample.

Panel (a) in (ref) shows how sensitive the natural rate estimates to the different starting dates are, particularly at the beginning of \cites{holston.etal:2017} sample, namely, from 1961 until about 1980, and at the end of the sample from 2009 onwards. Negative natural rate estimates are now also obtained for the U.S. when the beginning of the sample period is aligned with that of the Euro Area in 1972:Q1 (or 1967:Q1 which excludes the high GDP growth period). Since the natural rate is defined as the sum of trend growth $g_{t}$ and `other factor' $z_{t}$, we can separately examine the contribution of each of these factors to $ r_{t}^{\ast }$. From panel (b) in (ref) it is evident that the filtered trend growth estimates\ are the primary driver of the excessive sensitivity in $r_{t}^{\ast }$ over the 1961 to 1980 period. For instance, in 1961:Q1, these estimates can be as high as 6 percent, or as low as 3 percent, depending on the starting date of the sample. Also, the differences in these estimates stay sizeable until 1972:Q1, before converging to more comparable magnitudes from approximately 1981 onwards. Apart from the estimate using the very long sample beginning in 1947:Q1 (see the blue line in panel (b) of (ref)), the other four remain surprisingly similar, even during and after the financial crisis period, that is, from mid 2007 to the end of the sample in 2019:Q2. Thus, the `front-end' variability of the natural rate estimates are driven by the `front-end' variability in the estimates of trend growth $g_{t}$.

In panel (b) of (ref), I also superimpose MUE and MMLE (smoothed) estimates of trend growth from \cites{stock.watson:1998} model in (\oldref{eq:tvp_sim}), as well as (smoothed) estimates from the (correlated)\ UC model in (\oldref{Stage1:mod}) to provide long-sample benchmarks of trend growth from simple univariate models which can be compared to \cites{holston.etal:2017} estimates. These are the same estimates that are plotted in panels (b) and (c) of (ref). To avoid cluttering the plot with additional lines, I do not plot the mean and median estimates computed over the more recent expansion periods as was done in (ref). Note, however, that the MUE estimate overlaps with the mean and median of GDP\ growth from 2009:Q3 onwards and can thus be used as a representative for these model free `average' estimates of GDP growth since the end of the financial crisis. Comparing the Kalman Filter based estimates from the various starting dates to the MUE, MMLE, and UC (smoothed) ones shows how different these are, particularly, from 2009:Q3 until the end of the sample.\ In the immediate post-crisis period, the (one-sided) filter based estimates are `pulled down' excessively by the sharp decline in GDP and `converge' only slowly at the very end of the sample period towards the three long-sample benchmarks. Trend growth is severely underestimated from 2009:Q3 onwards, and this is reflected in the estimate of $r_{t}^{\ast }$.

In (ref) in the \hyperref[appendix]{Appendix}, I show plots of (real) GDP\ growth and the recursively estimated mean of GDP growth over the post financial crisis period from 2009:Q3 to 2019:Q2. Trend growth stays rather stable between $2\%$ and $3\%$ over nearly the entire period, settling at around $2.25\%$ in 2014:Q2 and remaining at that level. Moreover, it is never close to the filtered estimate of holston.etal:2017 from 2009:Q3 to 2014:Q3. In (ref) , I plot the mean as well as median 10 year ahead annual-average (real) GDP growth forecasts from the Survey of Professional Forecasters\ (SPF) from 1992 to 2020.\footnote{ The data were downloaded from: \url{https://www.philadelphiafed.org/research-and-data/real-time-center/survey-of-professional-forecasters/data-files/rgdp10} (accessed on the $27^{th}$ of July, 2020).} These forecasts also remain fairly stable between $2\%$ and $3\%$ from 2008 until 2017, and drift only marginally lower towards the very end of the sample. In (ref), Vanguard investor survey based 3 year and 10 year ahead expectations of (real) GDP growth from February 2017 to April 2020 are plotted. These are taken from Figure II on page 5 in giglio.etal:2020 . The 10 year expected growth rate shown in the right panel of (ref) fluctuates (mainly) between $2.8\%$ and $3.2\%$ (the 3 year expected growth rate on the left is somewhat lower). All three plots suggest that following the financial crisis, trend growth in GDP\ is unlikely to have dropped to the 1.3% estimate of holston.etal:2017.

Looking at the estimates of `other factor' $z_{t}$ in panel (c) of (ref), we can see that it is the end of the sample, namely, from 2009:Q1 to 2019:Q2, that is most strongly affected by the different starting dates.\footnote{ There is also some variability beginning in the 70s until the 80s, but this variation seems to be largely due to the noisier nature of the filtered estimates and is not visible from the more efficient smoothed estimates shown in panel\ (c) of (ref). The differentiation here is not important. The point to take away from this discussion is that the period following the financial crisis yields very different estimates from the two shorter samples, irrespective of whether smoothed or filtered estimates are used in the construction of the natural rate.} In particular the two $z_{t}$ estimates that are based on the shorter samples starting in 1967:Q1 and 1972:Q1, which exclude the `excessive' GDP growth period at the beginning of \cites{holston.etal:2017} sample, generate substantially more negative $z_{t}$ estimates. For instance, in 2009:Q1, the 1972:Q1 based $ z_{t}$ estimate is $-2.87$ while \cites{holston.etal:2017} is $-1.22$. Also, the $z_{t}$ estimates from the shorter samples are well below $-2$ over nearly the entire 2014:Q4 to 2019:Q2 period.\footnote{ This is even more pronounced in the smoothed estimates of $z_{t}$ shown in panel (c)\ of (ref).} What is particularly interesting to highlight here is how stable (and very close to zero) the estimates of $ z_{t} $ are from the four earlier sample starts from 1947:Q1 to about 1971:Q3. Given the rapid change in demographics and population growth, as well as factors related to savings and investment following the end of World War II, one would expect $z_{t}$ to capture this change. Even if we look at the period until 1990:Q1, apart from the noise in the estimates, no apparent upward or downward trend in $z_{t}$ is visible from panel (c) of (ref). Thus, the Baby Boomer generation entering the workforce shows no effect on $z_{t}$. Only from 1990:Q2 onwards is a decisive downward trend in the estimates of $z_{t}$ visible.

holston.etal:2017 initialize the state vector at zero for the $z_{t}$ elements of $\boldsymbol{\xi }_{t}$. This evidently has an anchoring effect on `other factor' $z_{t}$ at the beginning of the sample. In the model, it acts like a normalisation, as it implies that the natural rate is driven solely by trend growth $g_{t}$ at sample start. Although $z_{t}$ follows a (zero mean) random walk, so that an initialisation at zero seems appealing from an econometric perspective, this initialisation has an important impact on the economic interpretation of $z_{t}$ that should be more openly discussed if one is to view `other factor' $z_{t}$ as a factor relating to structural changes in an economy. Due to its large impact on the downward trend in the estimates of the natural rate, understanding exactly what $z_{t}$ captures and how the zero initialisation affects these estimates is crucial from a policy perspective.

One final point that needs to be raised relates to \cites{holston.etal:2017} preference for reporting filtered estimates of the latent states, as opposed to smoothed ones. It is well know that the mean squared error (MSE) of the filtered states will in general be larger than the MSE\ of the smoothed states (see the discussion on page 151 in harvey:1989). This is not surprising, as the smoothed estimates use the full sample --- and therefore more information --- to estimate the latent states, leading to more efficient estimates. Moreover, reporting filtered estimates `precludes' the use of a diffuse prior for the $I(1)$ state vector, since it generates extreme volatility in the filtered estimates of the states at the beginning of the sample period. This is not the case with the smoothed estimates. The large variability in the filtered states is particulary visible from the three quantities of interest, ie., the estimates of $ r_{t}^{\ast }$, $g_{t}$ and $z_{t}$, and less so from the output gap (cycle) estimates.

While it is frequently claimed that the filtered states are `real time' estimates, and are thus more relevant for policy analysis, one can see that this cannot be a valid argument in the given context. Not only are the parameter estimates of the model (ie., the estimates of $\boldsymbol{\theta } _{3}$ in (\oldref{eq:theta3})) based on full sample information, the GDP and PCE\ inflation data that go into the model are also not real time data, that is, data that were available to policy makers at time $t<T$. Reporting filtered (one-sided) estimates of the states as in holston.etal:2017 or as on the FRBNY website where updates are provided is undesirable from an estimator efficiency perspective.

\stepcounter{section} \Osubsection*{A.\arabic{section}. {Conclusion }} \addcontentsline{toc}{subsection}{\currentname}

\cites{holston.etal:2017} implementation of \cites{stock.watson:1998} Median Unbiased Estimation (MUE) in Stages 1 and 2 of their procedure to estimate the natural rate of interest from a larger structural model is unsound. I show algebraically that their procedure cannot recover the ratios of interest $\lambda _{g}=\sigma _{g}/\sigma _{y^{\ast }}$ and $\lambda _{z}=a_{r}\sigma _{z}/\sigma _{\tilde{y}}$ needed for the estimation of the full structural model of interest. \cites{holston.etal:2017} implementation of MUE\ in Stage 2 of their procedure is particularly problematic, because it is based on an `unnecessarily' misspecified model as well as an incorrect MUE\ procedure that spuriously amplifies their estimate of $ \lambda _{z}$. This has a direct and consequential effect on the severity of the downward trending behaviour of `other factor' $z_{t}$ and thereby the magnitude of the estimate of the natural rate.

Correcting their Stage 2 model and the implementation of MUE\ leads to a substantially smaller estimate of $\lambda _{z}$ of close to zero, and an elimination of the downward trending influence of `other factor' $ z_{t}$ on the natural rate of interest. The correction that is applied is quantitatively important. It shows that the estimate of $\lambda _{z}$ based on the correctly specified Stage 2 model is statistically highly insignificant. The resulting filtered estimates of $z_{t}$ are very close to zero for the entire sample period, highlighting the lack of evidence of `other factor' $z_{t}$ being important for the determination of the natural rate in this model. Obtaining an accurate estimate of trend growth for the measurement of the natural rate is therefore imperative. To provide other benchmark estimates of trend growth, I construct various simple alternative $g_{t}$ estimates and compare those to the estimate from holston.etal:2017. I find the latter one to be too small, particularly in the immediate aftermath of the global financial crisis.

Lastly, I\ discuss various other issues with \cites{holston.etal:2017} model that make it unsuitable for policy analysis. For instance, \cites{holston.etal:2017} estimates of the natural rate, trend growth, ` other factor' $z_{t}$ and the output gap are extremely sensitive to the starting date of the sample used to estimate the model. Using data beginning in 1972:Q1 (or 1967:Q1) leads to negative estimates of the natural rate as is the case for their Euro Area estimates. These negative estimates are again driven purely by the exaggerated downward trending behaviour of `other factor' $z_{t}$. The 1972:Q1 date was chosen to match the sample used in the estimation of the Euro Area model. Only the Euro Area estimates of the natural rate turn negative in 2013, and only the Euro Area sample starts in 1972:Q1 (the others start in 1961:Q1). The fact that it is possible to generate such negative estimates of the natural rate from \cites{holston.etal:2017} model for the U.S. as well by simply adjusting the start of the estimation period suggests that the model is far from robust, and therefore inappropriate for use in policy analysis.

Also, any Kalman Filtered (or Smoothed) estimates of the state vector will be a function of the observable variables that enter into the model. If the central bank controlled nominal interest rate is one of these observables, a confounding effect between $r_{t}^{\ast }$ and $i_{t}$ will arise, because any central bank induced change in the policy rate $i_{t}$ is mechanically transferred to the natural rate via the estimate of the state vector $ \skew{0}\boldsymbol{\hat{\xi}}_{t|t}$. This makes it impossible to answer `causal' questions regarding the relationship between $r_{t}^{\ast } $ and $i_{t}$, as one responds as a direct consequence to changes in the other.

{-5.4mm}

\ifthenelse{\equal{1}{1}}{ \IfFileExists{_figstabs.5e.tex}{\cleardoublepage \BAP\EAP}}

\ifthenelse{\equal{1}{1}}{ \cleardoublepage \changepage{0mm}{0mm}{0mm}{0mm}{0mm}{0mm}{0mm}{0mm}{5mm} \setcounter{page}{1} \setcounter{figure}{0} \setcounter{table}{0} \setcounter{equation}{0} \setcounter{footnote}{0} \setcounter{section}{0} \setstretch{1.234} \IfFileExists{_appendix.5e.tex}{\cleardoublepage \BAP\EAP}}

\cleardoublepage

\ifthenelse{\equal{1}{1}}{ \changepage{3mm}{0mm}{0mm}{0mm}{0mm}{-3mm}{0mm}{0mm}{0mm} \let\Osubsection\subsection \stepcounter{section} \Osubsection*{A.\arabic{section}. {R-Code Snippets }} \addcontentsline{toc}{subsection}{\currentname}

This sections shows various parts of the R-Code that is provided by holston.etal:2017 in the zip file from \url{https://www.newyorkfed.org/medialibrary/media/research/economists/williams/data/HLWCode.zip}. Below the code next to each of the headers, the name of the R-file is listed from which the code is displayed.

\multilstinputlisting{./HLW_Code_4tex/}{rstar.stage3.R}{1}{123}{R:stage3}

\cleardoublepage

\multilstinputlisting{./HLW_Code_4tex/}{calculate.covariance.R}{1}{31}{R:covar} \null

\cleardoublepage

\multilstinputlisting{./HLW_Code_4tex/}{unpack.parameters.stage1.R}{1}{45}{R:unpack} \null

\cleardoublepage

\multilstinputlisting{./HLW_Code_4tex/}{kalman.states.wrapper.R}{1}{45}{R:wrapper} \null

\cleardoublepage

\multilstinputlisting{./HLW_Code_4tex/}{rstar.stage2.R}{1}{116}{R:rstarS2} \null

\cleardoublepage

\multilstinputlisting{./HLW_Code_4tex/}{median.unbiased.estimator.stage2.R}{1}{60}{R:MUE2} \null

\cleardoublepage

}