Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
109,002 characters · 17 sections · 60 citation commands
Indirect Inference for Locally Stationary Models
Key words: semiparametric, locally stationary, indirect inference, state-space models
Journal of Economic Literature Classification:\ C13, C14, C22
Time-varying economic and financial variables, and relationships thereof, are stable features in applied econometrics. Notable examples include asset pricing models with time-varying features ghysels1998stable, wang2003asset and trending macroeconomic models stock1998median, phillips2001trending. While classical analyses of time series are built on the assumption of stationarity, data studied in finance and economics often exhibit nonstationary features.
Many different schools of modeling and estimation methods are used to accommodate the nonstationary behavior of observed time series data. In particular, statistical tools developed for locally stationary processes provide a convenient means of conducting analyses of trending economic and financial models. Heuristically, local stationarity implies that a process behaves in a stationary manner (at least) in the vicinity of a given time point but could be nonstationary over the entire time horizon. For certain widely-studied time series models, slowly time-varying parameters ensure local stationarity under some regularity conditions; for instance, see dahlhaus1996kullback and D97 (AR(1)), DSR06 (ARCH($\infty$)), dahlhaus2009empirical (MA($\infty$)), KL12 (Diffusion processes) and koo2015let (GARCH(1,1) with a time-varying unconditional variance) among many other classes of locally stationary processes.
While many classes of well-known time series models can be generalized to locally stationary processes, it is worth noting that estimation and inference procedures developed in one class of locally stationary processes often cannot be applied to a different class of locally stationary processes. In particular, many estimation methods for locally stationary processes are composed of estimation approaches that primarily focus on local regression with closed-form estimators, local maximum likelihood estimation (MLE) with a closed-form likelihood function (in the time domain) and spectral density approach (in the frequency domain), all of which could be intractable or simply difficult to implement for various locally stationary extensions of commonly used structural econometric models; we refer to vogt2012nonparametric, DSR06 and dahlhaus2009empirical, for examples. As such, model specifications compatible with the above statistical methods are rather limited and cannot be used for estimation and inference in more complicated locally stationary models, such as, for instance, models with latent variables or unobservable factors.
More importantly, structural models of economic and financial relationships commonly rely on the use of latent variables to represent information that is unavailable to the econometrician. This modeling approach implies, almost by definition, that simple (closed-form) representations for the conditional distributions of the endogenous variables are unavailable, with simple straightforward estimation methods often infeasible as a consequence. In such cases, if we were to extend common locally-stationary models to include the latent variables that are necessary to structurally model phenomena found in economics and finance, this would render the existing estimation methods used for such models infeasible. For instance, this situation arises in state-space models if either the measurement or state transition densities do not have closed forms, as in the case of stochastic volatility models. A secondary example is the fact that estimation of univariate locally stationary diffusion models cannot be straightforwardly extended to versions of these models with stochastic volatility.
To circumvent the above issue, and to help proliferate the use of locally stationary models and methods in econometrics and finance, we propose a novel nonparametric indirect inference (hereafter, II) method to estimate locally stationary processes. Instead of estimating complex structural locally stationary models directly, we indirectly obtain our estimator by targeting consistent estimators of simpler auxiliary models, and use these consistent estimates to conduct inference on the structural parameters. See, smith1993estimating, GMR93 and gourieroux1996simulation for discussion of indirect inference in parametric models.
To illustrate the main idea behind our nonparametric II approach for locally stationary processes, we consider the following motivating example. Suppose that the true data generating process evolves according to
where $\xi(t/T)>0,$ for all $t\leq T$. This locally stationary multiplicative stochastic volatility (LS-SV) model decomposes volatility into a short-term, latent volatility process, $h_t$, and a slowly time-varying component, captured by $\xi(\cdot)$, and can capture a wide range of volatility behaviors. The above model allows for non-stationary, but slowly changing, volatility dynamics, which may result from the transitory nature of the business cycle.
Suppose that we wish to estimate and conduct inference on the unknown volatility function $\xi(\cdot)$ in (ref). While (G)ARCH-based versions of the locally stationary volatility model have been analyzed by several researchers (see, e.g., DSR06, ER08, fryzlewicz2008normalized, and koo2015let), since the latent volatility process, $h_t$, pollutes the observed data, $Y_{t,T}$, it is not entirely clear how to estimate parameters in (ref). Indeed, largely due to this fact, locally stationary volatility models have not been previously explored in the literature, even though their stationary counterparts form the backbone of many empirical studies in finance and financial econometrics.
In this paper, we generalize the II approach of GMR93 to present a convenient estimator for unknown functions in locally stationary models, such as the LS-SV model. This approach to II estimation relies on a locally stationary auxiliary model that can be easily estimated using the observed data and that captures the underlying features of interest in the structural model. For example, in the context of the LS-SV model, a reasonable auxiliary model would be the locally stationary GARCH model:
where $\rho(t/T)>0$ for all $t\leq T$, and where $z_t$ is an error process.
The remainder of this paper further develops the ideas behind this estimation method in the context of a general locally stationary model and establishes the asymptotic properties of the proposed estimation procedure under regularity conditions. To establish the asymptotic properties of these II estimators, we must first develop conditions that guarantee locally stationary models admit consistent estimators of their corresponding limit values. This is itself a novel result since the vast majority of research into locally stationary models has focused on estimators defined by relatively simple criterion functions, and all under the auspices of correct model specification. Indeed, kristensen2019local is the only other study of which the authors are aware that treats genuinely misspecified locally stationary models. These new results for locally stationary estimators of the auxiliary model enable us to deduce the asymptotic properties of our proposed II estimator for the structural model parameters.
The estimation procedure proposed herein is demonstrated through two Monte Carlo examples, and an empirical application. The empirical application applies the LS-SV model to examine the volatility structure of several commonly analyzed Fama-French portfolios. We find that most of these portfolios display time-varying volatility patterns that broadly track the underlying (low-frequency) expansion and contractions of the United States economy.
The remainder of the paper is organized as follows. Section (ref) introduces the general model and the related framework. In Section (ref), we present our general approach and define the corresponding local II (L-II) estimators for a general locally stationary model. Section (ref) develops asymptotic results that demonstrate the properties of this estimation procedure. Simulation results for a simple example of a locally stationary moving average model of order one are discussed in Section (ref). In Section (ref) we analyze the locally stationary stochastic volatility model. We consider a small Monte Carlo to demonstrate our estimation method, then apply this method to analyze the volatility behavior of Fama-French portfolio returns, where we find ample evidence for smoothly time-varying nonlinear volatility dynamics over the sample period. All proofs are relegated to Appendix (ref). The tables and figures associated with the application in Section (ref) are given in Appendix (ref). The proof of Corollary 2 and additional details for the LS-SV model are provided in the supplementary appendix.
Throughout this paper, the following notations are used. The symbol $\mathbb{R}$ denotes the real numbers, while $\mathbb{N}$ denotes the natural numbers. For $x\in\mathbb{R}^d$, we let $\|x\|$ denote the Euclidean norm, while $|\cdot|$ denotes the absolute value function, and for $\Omega$ a $d\times d$ positive-definite matrix, we let $\Vert x\Vert^2_{\Omega}:=x'\Omega x$ denote the weighted norm of $x$. For $g:\mathbb{R}^{d}\rightarrow \mathbb{R}$ denoting a given function, we let $\|g\|_\infty:=\sup_{x\in\mathbb{R}^d}|g(x)|$ denote the sup-norm. For an unknown parameter $\theta$, the subscript $0$ denotes the true value of $\theta$. The quantities $O_p(\cdot)$ and $o_p(\cdot)$ denote the usual big $O$ and little $o$ in probability. $C$ denotes a generic constant that can take different values in different places.
We assume the researcher is interested in conducting inference on a model in the class of locally stationary processes.
The magnitude of $\eta$ captures the degree of approximation of $y_{t/T,t}$ to $Y_{t,T}$, which reflects the characteristics of the underlying processes of interest. The larger $\eta$, the better the approximation. We do not specify the magnitude of $\eta$ to maintain generality, which allows us to represent various types of processes, and instead allow $\eta$ to vary from model to model. See, for instance, DSR06 for ARCH($\infty$), KL12 for diffusion processes, vogt2012nonparametric for AR processes and dahlhaus2009empirical for MA processes among many other processes.
We consider that the process $\{Y_{t,T}\}$ is generated from the following {locally stationary} structural model:
where both $r(\cdot)$ and $\varphi(\cdot)$ are real-valued functions that are known up to the unknown function $\theta_0$. The function of interest is $\theta_{0}\in\mathcal{H}_\theta$, where $(\mathcal{H}_\theta,\|\cdot\|)$ denotes a normed vector space of function. The structural model, and $\theta_0$ satisfy the following regularity conditions.
The structural model in (ref) is quite general and can accommodate many interesting processes, including models with complex time-varying features, such as time-varying autoregressive conditional heteroskedasticity (ARCH). In addition, the structural model in (ref) can always be augmented with additional exogenous regressors at the cost of additional notation. Such regressors may be used, for instance, to capture some conditionally heteroskedastic features of the data. Critically for our purposes, under Assumption (ref), if $\theta_0(\cdot)$ were known, simulated realizations of $\{Y_{t,T}\}$ could easily be generated from the model in equation (ref).\footnote{We note here that Assumption (ref)(iii) is standard in the II literature. Indeed, GMR93 argue that this is not a real assumption since the error term “can always be considered as a function of a white noise with a known distribution and of a parameter which can be incorporated” into the unknown parameters.} If the process in (ref) is locally stationary, inference on $\theta_{0}(\cdot)$ can be carried out through an approximate structural model defining a stationary process indexed by $u\in\mathcal{U}$, where $\mathcal{U}$ denotes the domain of re-scaled time point $u=t/T$, i.e. $\mathcal{U}=[\delta,1-\delta]$ with a positive $\delta=o(1)$:
Lemma (ref) is consistent with Proposition 3.1 of dahlhaus2019towards, and implies that in the neighborhood of a re-scaled time point $u = t/T$, the local behavior of $\{Y_{t,T}\}$ can be approximated by the behavior of $\{{y}_{u,t}\}$. Consequently, statistical analysis on $\{Y_{t,T}\}$ can be based on a collection of locally stationary processes $\{y_{u,t}:u\in\mathcal{U}\}$.
Under local stationarity, we will demonstrate that estimation of the unknown (vector) function $\theta_{0}(\cdot)$ in (ref) can proceed through a local version of II (L-II) conducted at the time points $u=t/T$. This approach relies on the fact that, for any $u\in\mathcal{U}$, $\theta_0(u)$ in (ref) satisfies $\theta_{0}(u)\in\Theta\subset\mathbb{R}^{d_{\theta}}$; i.e., in the locally stationary structural model we view the function of interest as a map $\theta_{0}(\cdot):\mathcal{U}\mapsto\Theta$. The assumption that $\theta_{0}(\cdot)$ is our only parameter of interest is without loss of generality as we may always redefine $\theta_{0}(\cdot)$ to include those elements (time-varying or otherwise) of the distribution for the errors that are unknown. This paper is particularly concerned with estimation and inference when the structural model, (ref), rules out direct estimation approaches developed in the existing literature, for instance, due to the presence of latent variables that make computation of the likelihood function intractable.
Consider that our goal is to estimate the unknown map $\theta_0:\mathcal{U}\mapsto\Theta$ at a given point $u\in\mathcal{U}$. Since $\theta_0(u)\in\Theta\subset\mathbb{R}^{d_\theta}$, we associate to this unknown function (evaluated at the point $u$) a vector $\theta\in\Theta$. Even if the vector $\theta$ can not be estimated by direct means, since $\{y_{u,t}\}_{}$ is stationary (at the fixed value $u$) we can easily simulate a realization of this series by replacing $\theta_0(u)$ in equation (ref) by $\theta$. For fixed $u \in \mathcal{U}$ and some $\theta\in\Theta$, a simulated series $\{\tilde{y}_{u,t}(\theta)\}_{t\le T}$ can be generated according to
where $\tilde{\nu}_t$ denotes a simulated realization of the random variable $\nu_t$.\footnote{The use of slightly misspecified simulators in II is not uncommon, see, e.g., DGR2007, AEA2013, bruins2015, and frazier2018indirect for examples of misspecified simulators in the context of II estimation. In this sense, we follow the above papers in that the version of the structural model used to simulate data is a (locally) misspecified version of the true DGP.} Throughout the remainder, a tilde, $\tilde{}$, over a variable will denote that this variable is simulated and when no confusion will result we drop simulated series dependence on $\theta$, e.g., we take $\tilde{y}_{u,t}$ to mean $\tilde{y}_{u,t}(\theta)$.
Given the simulated series $\{\tilde{y}_{u,t}\}_{t\le T}$, II estimation of $\theta_0(u)$ can then proceed by minimizing the difference between statistics calculated from the observed data, $\{Y_{t,T}\}_{t\le T}$, and the simulated data, $\{\tilde{y}_{u,t}\}_{t\le T}$. Repeating this procedure at a collection of points $u_1,\dots,u_m$ would then yield an estimate of the unknown function $\theta_0(\cdot)$.
To employ our L-II estimation method, we specify an auxiliary model defined by the unknown (vector) function $\rho(\cdot)\in\mathcal{H}_\rho$, with $(\mathcal{H}_\rho,\|\cdot\|)$ a vector space of functions, and where $\rho(\cdot):\mathcal{U}\mapsto\Gamma\subset\mathbb{R}^{d_\rho}$ with $d_{\rho}\geq d_{\theta}$. Similar to the structural function of interest, for any given $u\in\mathcal{U}$ we associate to the unknown function $\rho(u)$ a vector $\rho\in\Gamma\subset\mathbb{R}^{d_{\rho}}$. In general, we will only emphasize the parameters' dependence on the point $u$ when necessary.
Reflecting the features of the true structural model, the auxiliary model is chosen such that it allows for direct estimation of $\rho(\cdot)$. We estimate $\rho(\cdot)$ at the point $u$, i.e., $\rho=\rho(u)$, by minimizing a local criterion function: for kernel function $K(\cdot)$ and bandwidth parameter $h$, define
where $g(\cdot)$ is a known function whose properties we later specify. Note that, technically $M_{T}[\rho;u]$ depends on the array $\{Y_{t,T}\}_{t\le T}$, however, we obviate this dependence to keep notation as simple as possible. Given $M_{T}[\rho;u]$, an estimator for $\rho(u)$ can be defined as
The explicit dependence of $\hat{\rho}(u;\theta_0(u))$ on $\theta_0(u)$ clarifies that the auxiliary estimator depends on the unknown $\theta_0(\cdot)$ at the point $u$. However, throughout the remainder, to simplify notation, we obviate this explicit dependence and simply define $\hat\rho(u):=\hat{\rho}(u;\theta_0(u))$.
It is natural to consider an auxiliary model which allows for simple estimation of the auxiliary parameters. One such useful class of auxiliary models will be nonlinear regression models of the type considered in robinson1991time and zhang2015time: for $Z_{t,T}$ a triangular array of variables that are measurable at time $t$, and exogenous with respect to the error term $\eta_t$, the auxiliary model is given as $$Y_{t,T}=f\left(Z_{t,T};\rho(t/T)\right)+\eta_t,$$ where $f(\cdot)\in\boldsymbol{F}$ is known, up to the unknown $\rho(\cdot)$, and where $$\boldsymbol{F}:=\{f:|f(x,\rho_1)-f(x,\rho_2)|\le b(x)\|\rho_1-\rho_2\|_{\infty},\;\rho_1,\rho_2\in \mathcal{H}_\rho\}.$$ The set $\boldsymbol{F}$ restricts the form of $f(\cdot)$ to be locally (in $x$) Lipschitz (in $\rho$), with this restriction being satisfied by many regression functions. Under this specific nonlinear regression model, $M_T[\cdot;u]$ could be the local least squares criterion
While nonlinear regression models are a useful class of auxiliary models, we do not wish to restrict our analysis solely to this class, and we therefore allow the criterion function $M_T[\rho;u]$ to be general. However, to ensure our theory can easily accommodate this case, we further specialize the structure of the auxiliary criterion function $M_{T}$: For some kernel function, $K(\cdot)$ and a bandwidth parameter, $h$, some known function $f(\cdot)\in\mathbf{F}$ and observable exogenous variables $Z_{t,T}$, we assume that
For $\{Y_{t,T}\}_{t\le T}$ denoting a set of observations from the locally stationary structural model (ref), satisfying Definition (ref), the auxiliary estimator $\hat\rho(u)$ in (ref) approximates the behavior of $\rho(\cdot)$ at the point $u$. Given $\hat\rho(u)$, an estimator of $\theta_0(u)$ can then be obtained by matching $\hat\rho(u)$ against a version that is calculated based on data simulated from the model under a given $\theta\in \Theta$, and a given $u\in\mathcal{U}$. However, we note that it is unclear in general how to simulate from the non-stationary structural model defined by (ref).
Therefore, instead of attempting to simulate from the model (ref), we invoke the local stationarity of $\{Y_{t,T}\}$ and generate (simulated) realization from the stationary process $\{y_{u,t}:u\in\mathcal{U}\}$, defined by (ref), which approximates $\{Y_{t,T}\}$ in the sense of Definition (ref). Such an II estimation approach is by construction “local” in that all we can recover is $\theta_{0}(u)$. An estimate of $\theta_0(\cdot)$ can be obtained by repeatedly applying this local II (L-II) approach at a given set of time points $\{u_{i}\}_{i=1}^m$, where $\max_{i}\Delta u_i=O(T^{-1})$ and $\Delta u_i:=u_{i}-u_{i-1}$.
More specifically, for some fixed $u_i\in\mathcal{U}$ and a corresponding candidate for $\theta_0(u_i)$, say, ${\theta}^{}={\theta}^{}(u_i)\in\Theta$, L-II then simulates data $\{\tilde{y}_{u_{i},t}^{}\}_{t\le T}$ from (ref) using simulated errors $\{\tilde{\nu}_t\}_{t\le T}$. Given $\{\tilde{y}_{u_i,t}^{}\}_{t\le T}$, we estimate the auxiliary parameters using
which corresponds to a simulated version of the local criterion function $M_T[\rho;u] $ in the vicinity of time point $u_i$. {Note that, similar to the notation we employ for $\hat{\rho}(u_i)$, the notation $\hat{\rho}(u_i;\theta^{})$ is an abbreviation for $\hat{\rho}(u_i;\theta^{}(u_i))$.}
Using $\hat{\rho}(u_i)$ and $\hat{\rho}(u_i;{\theta}^{})$, the L-II estimator of $\theta_{0}(u_i)$ can then be calculated, for positive-definite weighting matrix $\Omega$, as
Using the same simulated errors $\{\tilde{\nu}_{t}\}_{t=1}^{T}$, we may repeat the above procedure for $\{u_{i}\}_{i=1}^{m}$, with $0<u_{1}<u_{2}<\cdots<u_{m}<1$, {and $\max_{i}\Delta u_i=O(T^{-1})$}, to obtain an estimator of $\theta_{0}(\cdot)$.
The key feature of the above L-II procedure is that, due to the locally-stationary nature of (ref), the simulated series $\{\tilde{y}_{u_{i},t}\}_{t\le T}$ is stationary for each $u_{i}$, $i=1,...,m$. In this way, at each time point $u_{i}$, L-II matches a nonparametric estimator against a parametric estimator. As the following section illustrates, a consequence of this estimation approach is that the estimator $\hat{\theta}(\cdot)$ will inherit the asymptotic properties of the nonparametric estimator $\hat{\rho}(\cdot)$.
This section establishes the asymptotic properties of the L-II estimator. We establish the convergence (in probability) of $\hat{\theta}(\cdot)$ to $\theta_0(\cdot)$ and provide the asymptotic distribution of $\hat{\theta}(\cdot)$ under a fairly general setup.
Before presenting the details, we introduce the limit quantities that will be needed for our results. Consider the limit objective function and its minimizer corresponding to sample quantities, i.e. (ref) and (ref), such that, for $u\in\mathcal{U} =[\delta,1-\delta]$ and a small, positive $\delta=o(1)$, \[ \rho_0(u;\theta_0(u)):=\operatorname*{arg\,min}_{\rho\in \Gamma}\mathbb{M}_0[\rho;u],\text{ where } \mathbb{M}_0[\rho; u]:=\lim_{T\rightarrow \infty}E M_T[\rho;u]. \] When no confusion will result, we denote $\rho_0(u;\theta_0(u))$ by $\rho_0(u)$. The value $\rho_{0}(u)$ is the minimizer of the limit map $\rho \mapsto\mathbb{M}_0[\rho;u]$ and depends on the features of the true distribution and the true value of the unknown function, $\theta_0(\cdot)$, in the structural model.
Likewise, we require that the simulated auxiliary estimator has a well-defined probability limit. Recalling the stationary nature of the simulated data, $\tilde{y}_{u,t}$, such a requirement boils down to standard results for the consistency of quasi-maximum likelihood estimators for the pseudo-true value; see, e.g., white1982maximum and white1996estimation. The simulated counterpart to the pseudo-true parameter $\rho_0(u)$ is the map $\theta\mapsto\rho_0(u;\theta)$, which we define as
and where we remind the reader that we have suppressed the dependence of the simulated series $\tilde{y}_{u,t}$ on $\theta$ for notational simplicity.
To demonstrate the asymptotic properties of our proposed L-II approach, we employ the following regularity conditions.
For an arbitrary point $u\in\mathcal{U}$, define a local neighborhood of $\rho_{0}(u)$ as $\mathcal{E}:=\{\rho\in\Gamma:\|\rho-\rho_{0}(u)\|\le\varepsilon\}$ and $\mathcal{E}^{c}:=\{\rho\in\Gamma:\|\rho-\rho_{0}(u)\|>\varepsilon\}$.
Uniform (in $u$) consistency of the L-II estimator $\hat\theta(u)$ requires the uniform convergence of the auxiliary estimators $\hat{\rho}(u)$ and $\hat{\rho}(u;\theta)$ to their limit counterparts.
\color{black}
The (uniform) consistency of $\hat{\rho}(u)$ and $\hat{\rho}(u;\theta)$ allows us to deduce the uniform consistency of the L-II estimator.
In what follows, let $\Psi_{T}(\rho;u):=\sum_{t=1}^{T}q[Y_{t,T};f(Z_{t,T},\rho)]K\left(\frac{u-t/T}{h} \right)/Th$ and recall the definitions $\Psi_0(\rho;u):=(\partial/\partial\rho)\mathbb{M}_0[\rho;u]$ and $\mathcal{E}:=\{\rho\in\Gamma:\|\rho-\rho_{0}(u)\|<\varepsilon\}$. We deduce the asymptotic distribution of the L-II estimator under the following high-level regularity conditions.
In this section, we consider a simple generalization of the time-varying moving average model that allows the roots of the moving average lag polynomial to be time-varying. After presenting the model, we demonstrate how our L-II approach can be applied to estimate the model and present simulation results on the effectiveness of this strategy.
We consider the semiparametric locally stationary MA(1)-process
We further assume that $\epsilon_t$ is a white noise process with mean zero and unit variance, and $E\vert \epsilon_t\vert ^{4+\eta}<\infty$ for any arbitrarily small positive number $\eta$, and {we have that $\sup_{u\in \mathcal{U}}\vert\theta_0(u)\vert < 1$.}
Our goal is to estimate the unknown function $\theta_0(\cdot)$ {via our L-II approach.} In doing so, we approximate ((ref)) by a family of stationary MA(1) processes indexed by $u\in \mathcal{U}$ with some small trimming positive $\delta=o(1)$,
We consider an auxiliary model with the a locally stationary AR(1) structure:
For fixed $u$, the auxiliary model is a simple AR(1) model.
Recall that, in parametric MA models, when the roots lie near the region of non-inveribility, the resulting estimators can display a loss in accuracy. Therefore, since for any fixed $u$, the structural model is well-approximated by a parametric MA(1) model, it is likely that the same issue will be present if $\sup_u|\theta_0(u)|$ is close to unity.
We use the above auxiliary model to present a L-II estimator of $\theta_{0}(u)$. Algorithm (ref) describes the L-II estimation procedure for (ref).
In comparison with the general structure, the time-varying AR(1) auxiliary model in (ref) corresponds to taking $z_{u,t}=y_{u,t-1}$ and considering that $g(y_{u,t};\rho)=(y_{u,t}-\rho(u)y_{u,t-1})^2$. Note that it would also be possible to consider additional lags of $y_{u,t}$ in $z_{u,t}$ to accommodate LS-MA models of higher order. It is also useful to note that under weak conditions on the error term, the process $y_{t,T}$ defined in the auxiliary model (ref) is strong-mixing; see orbe2005nonparametric.
In this specific model, using the result of Corollary (ref), we can deduce the consistency result in Theorem (ref) to obtain the following uniform convergence of $\hat{\theta}(u) $ in the LS-MA(1) model to $\theta_0(u)$.\footnote{The proof of Corollary (ref) follows from Corollary (ref), however, for clarity we give a more primitive proof in the Supplementary appendix.}
We demonstrate the usefulness of the L-II approach using a series of Monte Carlo experiments. We consider a sample size of $T=1000$ generated according to the LS-MA(1) model $$Y_{t,T}=\epsilon_{t}+\epsilon_{t-1}\theta_{0}(t/T),\;\epsilon_{t}\sim \mathcal{N}(0,1).$$Data is generated according to one of three functional specifications for $\theta_{0}(u)$:
For inference on $\theta_{0}(\cdot)$, we use Algorithm (ref) with a Gaussian kernel and the rule of thumb bandwidth $h=1.06 T^{-1/5}$. We take $H=2$ for all simulation experiments.\footnote{As demonstrated in Theorem (ref), the choice of $H$ does not have an asymptotic impact on the estimates. However, in finite samples this choice may affect the estimated values of $\theta_0(\cdot)$, since a larger value of $H$ generally yields a smoother criterion function, and potentially a more accurate optimizer.} We estimate $\theta_{0}(\cdot)$ across the grid of points $u\in\{.05,.10,.20,\dots,.90,.95\}$.
We consider 5,000 replications of the above design across the three different specifications for $\theta_{0}(\cdot)$. The following three figures illustrate the sampling distribution, across the Monte Carlo replications for each of the three specifications.
Figure (ref) demonstrates the ability of the L-II approach to obtain consistent estimators of the unknown function $\theta_{0}(\cdot)$ over $u\in\{.05,.10,.20,\dots,.90,.95\}$ across the three Monte Carlo designs. The bounds are truncated due to the well-known boundary bias associated with local constant nonparametric estimation. We note that, outside of these bounds, given the relatively short nature of the time series, these estimators are likely to be poorly behaved. This issue can be addressed through the use of local-linear smoothing approaches.
The use of stochastic volatility to capture the conditional heteroskedastic movements of asset returns is now commonplace in economics and finance. Recently, however, several authors have suggested that volatility should be decomposed into short and long-run components (see, e.g, ER08 and engle2013stock). Such a decomposition has given rise to the class of multiplicative time-varying GARCH models, e.g. koo2015let. Such models decompose volatility into a short-run component, which is conveniently captured via a GARCH model, and a long-run component that slowly varies with larger macroeconomic factors that are captured nonparametrically.
The class of multiplicative GARCH models can capture both short and long-run features, however, it is generally accepted that stochastic volatility models are superior to GARCH models in terms of modeling flexibility and their overall ability to capture fluctuations in short-run volatility. Given this feature, one would suspect that a multiplicative extension of the standard SV model should perform well in many cases. While such a model would be similar to multiplicative GARCH models, the introduction of latent stochastic volatility ensures that direct estimation approaches become infeasible. However, this issue is immaterial for our L-II estimation approach since we can simulate the latent volatility
To this end, in this section we propose a new model where volatility evolves as the product of a short and long-run component: the long-run component is captured by a slowly time-varying function, and the short-run component is captured via an autoregressive SV model. In the context of simulation experiments, we demonstrate that our L-II approach can accurately estimate this new model. We then apply this model to analyze the volatility of monthly returns on twenty-five Fama-French portfolios, with the results indicating that long-run volatility changes dramatically over the sample period under analysis.
Given the general nature of this paper, we leave a thorough discussion on the theoretical properties of this new SV model for future study.
We now consider a multiplicative extension of the traditional stochastic volatility model. The observed demeaned data is generated according to
and where
with $\gamma_{\nu}$ the correlation coefficient between $\nu_{1,t}$ and $\nu_{2,t}$. In this model, the long-run trend is captured by the deterministic function $\sqrt{\xi({t/T})}$ whereas the short-run dynamics, $h_t$, are represented by the stochastic volatility model. We implicitly assume that $\{Y_{t,T}\}$ changes smoothly over time and if it were not for $\xi(\cdot)$, the slowly time-varying long-run trend, then $\{Y_{t,T}\}$ would be stationary. That is, we implicitly maintain that $\xi(\cdot)$ is uniformly positive and twice continuously differentiable, and $h_t$ is stationary, so that the process $\{Y_{t,T}/\sqrt{\xi(t/T)}\}_{}$ would be stationary. In the supplementary material, we give precise conditions on the function $\xi(\cdot)$ and the remaining parameters that ensure the resulting model is locally stationary.
Directly estimating the structural model (ref), and conducting statistical inference on the resulting estimates, is generally infeasible with existing methods. Instead, we propose to conduct inference on the structural model through L-II and by using as our auxiliary model the following locally stationary multiplicative GJR-GARCH model:
where $z_t \overset{iid}{\sim}\mathcal{N}(0,1)$ and $I_t = 0$ if ${y_{u,t}}/{\sqrt{\tau(u)}} \geq 0$, and $I_t = 1$ if ${y_{u,t}}/{\sqrt{\tau(u)}}< 0$. In this setting, we will use the parameters in the auxiliary model, $\rho(\cdot) = (\tau(\cdot),\omega,\alpha,\beta,\gamma)'$, to conduct inference on the parameters of interest in the structural model, $\theta(\cdot) = (\xi(\cdot),\mu,\phi,\gamma_{\nu},\sigma)'$.
koo2015let demonstrate that locally stationary multiplicative GARCH models can be estimated relatively easily. Note, however, that the symmetry of a GARCH(1,1) model would ensure that it is an unsuitable auxiliary model, as there is no parameter that can be readily matched to the correlation coefficient $\gamma_{\nu}$. Therefore, we employ the GJR-GARCH(1,1) model so that the leverage effect $\gamma_{\nu}$ is captured by the asymmetry parameter $\gamma$ in the auxiliary model.
Before we discuss estimation of the LS-SV model, we note that, due to the multiplicative nature of the model for $Y_{t,T}$ in (ref), an additional identification restriction is required in order to identify the unknown parameters. The restriction can be imposed on either the long-run or the short-run part. For instance, while koo2015let impose a restriction on the long-run component, engle2013stock impose a restriction on the short-run component. For our L-II, we impose a restriction on the short-run component for the LS-SV model because the L-II is applied over a finite number of fixed time points and therefore, a restriction on the long-run component in the structural model is difficult to implement.
In particular, we impose the restriction that $\mu = 0$ for the structural model. Equivalently, for the auxiliary multiplicative GJR-GARCH model, we restrict $\omega = 1-\alpha-\beta -\frac{\gamma}{2}$ such that the GJR-GARCH process has unit unconditional variance ($\frac{\omega}{1-\alpha-\beta-\frac{\gamma}{2}}$ = 1). Under this setup, we conduct our L-II as follows.
\paragraph{Estimation of the auxiliary model:} Using the observations $\{Y_{t,T}\}_{t\le T}$, we estimate the auxiliary multiplicative GJR-GARCH model \`{a} la ER08 and koo2015let. Specifically, from (ref), for $\mathcal{I}_t$ denoting the information set at $t$,
under the stationarity of $\sigma^2_t$ and $z_t$ and $\tau^{\ast}(u)=\tau(u)\exp(C)$ with $C=E(\log \sigma^2_t z_t^2|\mathcal{I}_{t-1})$.
We obtain an initial estimate $\log\hat{\tau}^{\ast}(u)$ as \[ \log \hat{\tau}^{\ast}(u) = \operatorname{arg min}_{\tau^{\ast} \in \mathbb{R}_{+}}\sum_{t=1}^{T}(\log y^2_{u,t}-\log \tau^{\ast}(u))^2 K_h(u-t/T), \] where $K_h(\cdot)=K(\cdot/h)/h$ with a bandwidth $h$. Once we obtain $\hat{\tau}^{\ast}(u)$, we calculate the intermediate estimator $\check{\tau}(u)$: \[ \check{\tau}(u)=\frac{\hat{\tau}^{\ast}(u) }{\int_{0}^{1}\hat{\tau}^{\ast}(u) du} \] because \[ \frac{{\tau}^{\ast}(u) }{\int_{0}^{1}{\tau}^{\ast}(u) du} = \frac{{\tau}(u)\exp(C)}{\int_{0}^{1}{\tau}(u)\exp(C) du}=\tau(u), \] when we impose a restriction that $\int_0^1\tau(u)du=1$.
{Note that the restriction, $\int_0^1\tau(u)du=1$ is not a model restriction but rather an estimation restriction that can be re-normalized or reconstructed arbitrarily. Once $\check{\tau}(u)$ is obtained, we estimate the GJR-GARCH parameters via maximum likelihood estimation based on the following transformed data $\check{y}_{u,t}=y_{u,t}\big{/}\sqrt{\check{\tau}(u)}$ and obtain the estimators $(\check{\omega},\check{\alpha},\check{\beta},\check{\gamma})'$. However, note that $\check{\rho}=(\check{\tau}(\cdot),\check{\omega},\check{\alpha},\check{\beta},\check{\gamma})'$ does not satisfy the restriction $\omega = 1-\alpha-\beta -\frac{\gamma}{2}$. To obtain a vector of parameter estimates that satisfy this restriction, we calculate $\hat{\tau}(u)=\check{\tau}(u)\left({\check{\omega}}/{1-\check{\alpha}-\check{\beta}-\frac{\check{\gamma}}{2}}\right)$ and use $\hat{\tau}(u)$ to construct $\hat{y}_{u,t} = y_{u,t}\big{/}\sqrt{\hat{\tau}(u)}$. Estimating the parameters in the GJR-GARCH model using the transformed dataset $\{\hat{y}_{u,t}\}_{t\le T}$ then yields $(\hat{\omega},\hat{\alpha},\hat{\beta},\hat{\gamma})'$. The vector of estimates $\hat{\rho}=(\hat{\tau}(u),\hat{\omega},\hat{\alpha},\hat{\beta},\hat{\gamma})'$ is then used in L-II as the auxiliary parameter estimates.\footnote{Imposing a restriction in maximum likelihood estimation is usually difficult but we avoid complicated constrained optimization in this way. This restriction or constraint is important for the L-II of this particular model. Another type of constraint is required for another type of structural and auxiliary models for L-II. We believe that imposing a general type of constraint in the context of L-II will open up another important research topic. We leave the analysis of constrained L-II for future research.}
\paragraph{Simulation of the structural model:} Based on (ref), for a given $u\in\mathcal{U}$, we simulate $H$ independent structural processes under the restriction $\mu=0$, for some value of $\theta\in\Theta$ according to:
with
In the simulation step, we restrict $\mu = 0$ to impose unit unconditional variance for the multiplicative SV model, which is compatible with the restriction on the auxiliary model, $\omega = 1-\alpha-\beta -\frac{\gamma}{2}$.
\paragraph{Estimation of the simulated structural model via the auxiliary model and L-II:} {For a given time point $u\in\mathcal{U}$, based on the simulated data $\{\tilde{y}^{[j]}_{u,t};j=1,...,H\}$, we first obtain a set of estimators $\{\hat{\rho}^{[j]}(u;\theta)\}_{j=1}^{H}$. Note that when $\{\hat{\rho}^{[j]}(u;\theta)\}_{j=1}^{H}$ is estimated for each fixed time point, $u$, the parameter $\tau(u)$ in the auxiliary model is an unknown constant, not a function. This implies that we just estimate the GJR-GARCH model based on the simulated data $\{\tilde{y}^{[j]}_{u,t};j=1,...,H\}$, to obtain $\{\check{\omega}^{[j]},\check{\alpha}^{[j]},\check{\beta}^{[j]},\check{\gamma}^{[j]}\}_{j=1}^{H}$ and then obtain $\{\hat{\tau}^{[j]}(u)\}_{j=1}^{H}$, such that $\hat{\tau}^{[j]}(u)=\frac{\check{\omega}}{1-\check{\alpha}-\check{\beta}-\check{\gamma}/2}$ thanks to the restriction $\omega = 1-\alpha-\beta -\frac{\gamma}{2}$. Then we create transformed or normalized data $\hat{y}_{u,t}=\tilde{y}^{[j]}_{u,t}\big{/}\sqrt{\hat{\tau}^{[j]}(u)}$ and obtain $\{\hat{\omega}^{[j]},\hat{\alpha}^{[j]},\hat{\beta}^{[j]},\tilde{\gamma}^{[j]}\}_{j=1}^{H}$. From $\{\hat{\rho}^{[j]}(u;\theta)\}_{j=1}^{H}$ we can then construct $\hat{\rho}(u;\theta)=\sum_{j=1}^{H}\hat{\rho}^{[j]}(u;\theta)/{H}$.}
{Based on $\hat{\rho}(u)$ and $\hat{\rho}(u;\theta)$, we search for the best candidate for the given time point $u$ and define the estimator $\hat{\theta}(u)$ as the solution to: $\operatorname*{arg\,max}_{\theta\in \Theta}-\|\hat{\rho}(u)-\hat{\rho}(u;\theta)\|^2_{\Omega}$ where $\Theta$ is the parameter space for $\theta_0(u)$. The above procedure can then be repeated across a grid of points, say $\{u_i\}_{i=1,...,m}$ to estimate the whole functional form of $\theta_0(\cdot)$. }
Summing up, Algorithm (ref) is employed for the L-II estimation of the locally stationary multiplicative stochastic volatility model.
We now conduct a Monte Carlo experiment to illustrate L-II estimation of the locally stationary multiplicative stochastic volatility (LS-SV) model . We fix the sample size to be $T=200$, and we generate 5000 Monte Carlo replications from the LS-SV model in equation (ref) with parameters values given by $$\mu=0,\phi=0.2,\gamma_{\nu}=-0.5,\sigma=1,$$ and where the long-run volatility component is given by $$\xi(t/T)=0.2\sin(0.5\pi t/T)+0.8\cos(0.5\pi t/T).$$ We take as our auxiliary model for this Monte Carlo experiment the LS-GJR-GARCH(1,1) auxiliary model in equation (ref).
Similar to the Monte Carlo experiments for the LS-MA(1) model, we estimate the auxiliary parameter via local constant estimation with a Gaussian kernel and rule of thumb bandwidth. We again set the number of simulations to be $H=2$. For full details of the estimation procedure, please refer to Algorithm (ref). Across each Monte Carlo replication we apply the LS-II approach, and record the estimated function $\hat{\xi}(\cdot)$.\footnote{Results for the parametric components of the model are similar to those obtained for other II estimators, and are not presented for the sake of brevity.} The estimation results for the unknown function are presented graphically in Figure (ref). Similar to the results for the LS-MA(1) model, the LS-II procedure yields good estimates of the unknown function.\footnote{Similar to the previous Monte Carlo, we truncate the function estimate due to boundary bias problems associated with the local-constant smoothing approach considered in this implementation.}
Herein, we analyse the behavior of monthly returns from January 1952 until December 2018 on 25 Fama-French portfolios formed from the intersection of five portfolios on size and five portfolios on book-to-market, and where the breakpoints for the portfolios are taken from the NYSE quintiles and are ordered from smallest to largest.\footnote{The data is freely available from Kenneth French's website.} The monthly return series on the Fama-French portfolios covers a long period of observation, and it is unlikely that these series display constant conditional covariance features over the entire sample period. In particular, while it is fairly widely accepted that these portfolios seem to display constant mean dynamics, the large fluctuations in the volatility of these series do not engender confidence that the conditional variance is constant throughout the sample period.\footnote{Considering an ARCH test of the demeaned returns for each of the 25 portfolios, where each test uses five lags, returns overwhelming support for the alternative hypothesis across all portfolios. The specific values can be found in the supplementary appendix.}
Moreover, given the long time-span over which the data is measured, we argue that it is not realistic to assume that the volatility dynamics that were present in the 1950s have persisted unchanged until 2018. In particular, it is likely that underlying macroeconomic factors would cause these portfolios to exhibit patterns of volatility that display both short-term and long-run fluctuations, which can not be adequately captured by a stationary volatility model. To capture the long-run volatility patterns in the data, we consider a LS-SV version of the Fama-French three factor model. For $r_{t,j},\;j\in\{1,...,25\}$, denoting excess returns on the $j$-th portfolio, we assume that $r_{t,j}$ evolves according to
where $r_{t,m}$ denotes excess returns on the market factor, $\text{SMB}_t$ is the size factor, and $\text{HML}_t$ is the value factor. We model the short-term volatility component $h_{t,j}$ as
where we require that the mean of the short-term SV component be zero to ensure the scale of $\xi(\cdot)$ can be properly identified. The above LS-SV model considers that volatility is the composition of two components: a long-run volatility trend that moves slowly and is captured by $\xi_j(t/T)$, and a term, measured by $h_{t,j}$, that captures short-term fluctuations around $\xi_j(t/T)$.
Estimation in the above LS-SV model can be carried out in two steps: first, we estimate the regression parameters to obtain $\hat{\alpha},\hat{\beta}_1,\hat{\beta}_2,\hat{\beta}_3$; in the second step, the residuals $$y_{t,j}=\left(r_{t,j}-\hat{\alpha}-\hat{\beta}_{1}r_{t,m} -\hat{\beta}_{2}\text{SMB}_{t}-\hat{\beta}_{3}\text{HML}_{t}\right)$$ are used within the L-II algorithm for the LS-SV model, along with a LS-GJR-GARCH auxiliary model (we refer the reader to Algorithm (ref) for specific implementation details). Before moving on, we note that the two-step nature of the L-II approach in this example means that it is straightforward to treat more complicated regression functions, such as, for instance, models with time-varying $\alpha$ and $\beta$. We refer the interested reader to the supplementary appendix where we consider an alternative specification for the conditional mean function that allows $\alpha,\;\beta$ to be time-varying.\footnote{These results largely mirror those given in the main text, and so we relegate these details to the supplementary material. In particular, we find that time varying versions of $\alpha$ and $\beta$ do not meaningfully deviate from constants for the sample period under analysis.}
L-II is used to estimate the short-term and long-run volatility components for all 25 portfolios. However, given the nature of the above estimation approach, uncertainty quantification is carried out using the local block-bootstrap (LBB) of paparoditis2002local. The LBB is operationally similar to the block bootstrap but accounts for the changing stochastic structure of the observation process. Given observed data $y_{1},\dots,y_{T}$ the LBB generates a bootstrapped series of data, $y_{1}^*,\dots,y_{T}^*$, via the following steps.
In the following examples, across each of the 25 portfolios, we implement the LBB using $R=999$ bootstrap replications. Furthermore, we set the LBB block size, $b$, to be $b=10$, and take the local bootstrap parameter, $B$, to be $B\approx0.11$.
The estimation results for $\alpha$ and $\beta$ are given in Table (ref), and the results for the parametric SV components are given in Table (ref). Focusing on the values of $\alpha,\;\beta$, we see that these estimated parameters are generally statistically significant and have the anticipated signs. Analysing Table (ref), we see that the short-term volatility parameters generally have statistically significant autocorrelation coefficients between 0.5 and 0.7, which indicates a moderate amount of short-term volatility persistence. The majority of the estimated values for $\sigma_v$ are between 1.5 and 2.0, indicating a relatively large level of noise in the short-term volatility process. Interestingly, none of the estimated leverage effects are statistically significant for the short-term volatility process. To ensure that this insignificance is not an artifact of the chosen auxiliary model, in Table (ref) we report 99% confidence intervals for the corresponding LS-GJR-GARCH auxiliary parameter $\gamma$, which captures the impact of asymmetric news on volatility, and where the confidence intervals are calculated using QMLE sandwich form standard errors. For 24 out of the 25 portfolios, the resulting LS-GJR-GARCH asymmetry parameter is statistically insignificant at the one percent significance level.
Leverage effects account for asymmetric reactions to volatility, possible due to larger macroeconomic forces. By their very nature, these macroeconomic forces are generally slowly varying, and their impact on volatility can then be adequately captured using the time-varying volatility approach considered herein. The insignificance of the estimated leverage effects can then be interpreted as follows: by decomposing volatility into a short-term and long-run component, and by modeling the impact of such macroeconomic forces nonparametrically, the leverage effect is soaked-up by the long-run volatility component; its inclusion in the short-term volatility component is then redundant and, hence, statistically insignificant.
We present the estimates of $\xi(\cdot)$ graphically in Figures (ref)-(ref) in the appendix. The reported confidence bounds are the corresponding pointwise, for each value of $u=t/T$, confidence bounds obtained using the LBB.
The long-run volatility component captures gradual changes in volatility, possibly due to slowly-varying macroeconomic factors that affect returns (see, e.g, ER08 and engle2013stock for a detailed discussion). Given this aim, the results in Figures (ref)-(ref) are compelling as they closely align with the larger macroeconomic risk profile of returns over the sample period under analysis. In particular, during the 1950s to the early 1960s most series display relatively low volatility that is either flat or slightly increasing till the early-to-mid 1960s, with the overall trend of most series decreasing after about 1965. This overall trend is then maintained all the way through the great moderation of the 1980s. However, after the end of the great moderation, virtually every series exhibits a significant upswing in long-run volatility. This pattern then continues and culminates around the time of the global financial crisis in the late 2000s, after which there is another sustained decrease in long-run volatility.
Given how well our results correspond to the overarching long-run volatility patters, we note that more than half of these return series now exhibit an additional steeping of long-run volatility. This may indicate that since 2016 we have entered into a new period of long-run macroeconomic volatility.
We propose a novel indirect inference estimator for locally stationary processes and thereby extend, for the first time, the use of indirect inference estimation to general classes of semiparametric models with slowly time-varying parameters. As part of this study, we also propose a novel local stationary multiplicative stochastic volatility (LS-SV) model. We leave two important topics for future research: the efficiency of the L-II estimator, and the ensuing semiparametric efficiency bound for the class of locally stationary models considered in this paper; and the incorporation of shape restriction for nonparametric estimation within L-II, which may improve efficiency, e.g. horowitz2017nonparametric, at the cost of a more complicated estimation approach.