Fitting Dynamically Misspecified Models: An Optimal Transportation Approach

refsection\begin{abstract} This paper considers filtering, parameter estimation, and testing for potentially dynamically misspecified state-space models. When dynamics are misspecified, filtered values of state variables often do not satisfy model restrictions, making them hard to interpret, and parameter estimates may fail to characterize the dynamics of filtered variables. To address this, a sequential optimal transportation approach is used to generate a model-consistent sample by mapping observations from a flexible reduced-form to the structural conditional distribution iteratively. Filtered series from the generated sample are model-consistent. Specializing to linear processes, a closed-form Optimal Transport Filtering algorithm is derived. Minimizing the discrepancy between generated and actual observations defines an Optimal Transport Estimator. Its large sample properties are derived. A specification test determines if the model can reproduce the sample path, or if the discrepancy is statistically significant. Empirical applications to DSGE models, affine term structure models, and trend-cycle decomposition illustrate the methodology and the results. \end{abstract} JEL Classification: C11, C12, C13, C32, C36.\newline Keywords: Semiparametric estimation, Model Evaluation. \baselineskip=18.0pt \thispagestyle{empty} \setcounter{page}{0} \section{Introduction} Structural estimation is routinely used to evaluate Economic theories and conduct counterfactual analyses with non-experimental data. As noted by domowitz1982, to make the analysis tractable - either analytically or numerically - some simplification is needed so that the model merely approximates the actual, more complex, data-generating process. Though misspecified, the model still provides tractable insights about causal mechanisms that can be used for policy evaluation and conduct counterfactual experiments. This paper is specifically interested in multivariate models of the form: \begin{align} y_t &= g(z_t,v_t;\theta), \quad z_t = h(z_{t-1},v_t;\theta), \end{align} where $y_t$ are observed variables such as output or inflation, $z_t$ are unobserved variables such as productivity, and $v_t$ are structural innovations. The functions $g$ and $h$ are known, or solved numerically, up to parameters of interest $\theta$. This is known as a state-space, or hidden Markov model. Examples include DSGE and structural asset pricing models. It is common to fit the model using a filtering algorithm that recovers the latent variables $z_t$ -- the Kalman or particle filter -- and then maximize the likelihood, or sample Bayesian posterior draws. This paper shows that several issues can arise when the dynamics in $g$ or $h$ are misspecified. First, for a given value $\theta$, filtered variables may not satisfy the model restrictions given by $g$ and $h$. For instance, the same average of a filtered series can be substantially non-zero even though the model describes a mean-zero process. This can make inferences on policy-relevant variables difficult, e.g., output gap or natural rate of interest, since their interpretation is model-dependent. Similarly, filtered structural shock processes can be cross-correlated even though the model specifies them as independent, a well-documented phenomenon in the DSGE literature. Second, likelihood estimates $\hat{\theta}_n$ might be hard to interpret since these coefficients may not characterize the dynamics of the filtered variables: the serial correlation of a Kalman filtered shock can differ substantially from its model-implied value. These two points are illustrated using a medium-scale DSGE model. Third, the likelihood is not defined when there are fewer structural shocks than observables. This limits the potential for model-based dimension reduction where a few structural shocks are used to summarize the comovement of financial or economic variables. This is also illustrated in the applications. This paper considers an optimal transportation approach to fitting dynamics models in ((ref)) with a mean-squared criterion. Fitting here refers to filtering latent variables and estimating the parameters of interest. The basic idea is to first flexibly approximate the true dynamics of the data; using a reduced-form model. Then, for a given $\theta$, a new sample is recursively constructed for which ((ref)) holds. At each time iteration, the procedure maps the observations to model-consistent data by transporting from the reduced-form to the model-predicted conditional distribution, i.e. the conditional mean-squared error between the sample and the model-consistent data is minimized. Finding the least-squares difference between the original and the new sample produces optimal transport estimates for the parameters of interest. A by-product is an optimal transport filtered series for $z_t$. The new data, filtered values, and estimates preserve model dynamics by construction, and are thus internally valid. Note that, unlike with i.i.d. data, the dependence structure here requires a different approach to implementing the transport. This is reflected in the use of a flexible reduced-form model and the iterative nature use of a conditional transport, where each step builds on the previous ones and is performed as many times as the sample size. Although there has been much progress in the computation of non-linear filters and numerical optimal transportation, the generic approach described above can be computationally demanding for larger models. Specializing to linear processes, a plugin rule for the optimal transport map is derived leading to closed-form expressions. The true dynamics of the data are approximated using a sieve approach through a vector autoregression of increasing order. The resulting algorithm has closed-form, is easy to implement, and numerically inexpensive. The associated estimator is semiparametric, as only the linear dynamics are specified. The closed-form plugin map extends to a class of non-linear models, though not as general as ((ref)). Because the transport map reduces to the identity map under correct model specification, the framework encompasses correctly specified structural models as a special case. For stationary linear processes, we derive the large sample frequentist results for the optimal transport estimator $\hat{\theta}_n$ which minimizes the mean-squared difference between original and model-consistent samples. Under standard regularity conditions, the estimates are shown to be consistent and asymptotically normal at a $\sqrt{n}$-rate. An expression for the asymptotic variance is derived under correct specification and misspecification. A specification test based on the mean-squared discrepancy between the sample paths is proposed and studied. The method and results cover a large class of models with an infinite moving average representation which includes linear state-space models. While the results are confined to frequentist estimation, it can also be of interest to extend the framework and consider quasi-Bayesian posterior sampling and inference. This goes beyond the scope of this paper. In Machine Learning, the sample Wasserstein distance between distributions is a popular tool for data analyses by optimal transportation. It is generally intractable, non-smooth, suffers from a curse of dimensionality, and is computationally demanding for estimation. This can limit its appeal for estimating models with a moderate or large number of observables and parameters. In the scalar case, the minimum Wasserstein distance estimator is fully parametric but has non-standard limiting distribution, which complicates inference. Entropic regularization is a popular way to circumvent some limitations of the Wasserstein distance, but it introduces bias. In contrast, the setting here is semiparametric -- only first and second order moments of the data and model are involved. The auxiliary model provides these moments for the data. We show that this allows for a computationally trivial closed-form solution to the transport problem, even for medium-scale DSGE models. The closed-form map is smooth, making estimation regular with $\sqrt{n}$-asymptotically Gaussian estimates. Three empirical applications illustrate the methodology and the issues discussed above on well-known models. First, a small-scale DSGE model from lubik2004 is estimated. The fit for inflation is rejected at the 5% level, which corroborates previous findings. Second, the medium-scale smets2007 model is estimated. Unlike previous studies, the fit for consumption is rejected as it does not match volatility and persistence. Further, the Kalman filtered series display irregularities consistent with misspecification, exacerbated by persistence in the data. Next, an affine term structure model based on hamilton2012 illustrates the dimension-reduction aspect: 3 factors explain most of the variation of 6 yields ranging from 1 month to 5 years. Finally, a trend-cycle decomposition shows that the optimal transport filter produces an interpretable cycle, while the Kalman filter produces a cycle estimate that is systematically positive over 1965-2008. \section{Filtering under Dynamic Misspecification} This section motivates the proposed Optimal Transport Filter (OTF). It illustrates the Kalman Filter (KF) under misspecification using a small DSGE model, then discusses the effects of misspecification on filtering in a general setting. \subsection{Motivating Example: A Small DSGE Model} The following illustrates two issues with Kalman filtering under misspecification: (i) filtered variables need not satisfy model constraints, (ii) the KF fits observables perfectly. In constrast, the OTF introduced in this paper: (i) recovers latent variables that satisfy model constraints, and (ii) indicates lack of fit in observables, which is informative about the misspecification. The model is taken from lubik2004: \begin{align} \begin{split} &y_{t} =E_{t}y_{t+1}-\tau (r_{t}-E_{t}\pi _{t+1})+g_{t}, \quad \pi _{t} =\beta E_{t}\pi _{t+1}+\kappa (y_{t}-z_{t}), \quad g_{t} =\rho _{g}g_{t-1}+\varepsilon_{gt}, \\ &z_{t} =\rho_{z}z_{t-1}+\varepsilon _{zt}, \quad r_{t} = \rho _{r}r_{t-1}+(1-\rho _{r})\psi _{1}\pi _{t} +(1-\rho _{r})\psi _{2}(y_{t}-z_{t}) +\varepsilon _{rt} \end{split}, \end{align} where $y_{t},\pi _{t}$, and $r_{t}$ are log-deviations of output, inflation, and the nominal interest rate from their steady states. The shocks $\varepsilon _{rt},$ $\varepsilon _{gt},\varepsilon _{gt}$ are iid Gaussian with mean zero and variances $\sigma_{r}^{2},\sigma _{g}^{2},\sigma _{z}^{2}$; and $\varepsilon _{gt}$ and $\varepsilon _{zt}$ are correlated with correlation $\rho_{gz}$. Using empirical estimates from Section (ref) below, $n = 10^6$ observations $(\tilde{y}_t,\tilde{\pi}_t,\tilde{r}_t)$ are simulated from model ((ref)). The KF and OTF are applied to the data, without re-estimating the model to focus on filtering. \\ Correct specification. As a benchmark, the KF and OTF are applied to the correct specification ((ref)) using the true parameter values. The $R^2$ between the filtered observables $({y}_t,{\pi}_t,{r}_t)$ and the data $(\tilde{y}_t,\tilde{\pi}_t,\tilde{r}_t)$ equals one for the KF and $0.9999$, or higher, for the OTF.\footnote{$R^2 = 1-\text{ESS}/\text{TSS}=1- \sum_{t=1}^n (\tilde{y}_t-y_t)^2/\sum_{t=1}^n (\tilde{y}_t-\tilde{y})^2$ with $\tilde{y}$ the sample mean of $\tilde{y}_t$. } The three $R^2$ are similar when comparing the true $(\tilde{\varepsilon}_{rt},\tilde{g}_{t},\tilde{z}_{t})$ with the filtered $(\varepsilon_{rt},g_{t},z_{t})$ shocks. As expected, when the model is correctly specified, both filters recover the latent variable accurately and show a perfect fit for the observables.\\ Misspecification. Now suppose a researcher misspecifies the model by assuming the central bank targets expected, rather than current, inflation: \begin{align} r_{t} = \rho _{r}r_{t-1}+(1-\rho _{r})\psi _{1}\mathbb{E}_t(\pi _{t+1}) +(1-\rho _{r})\psi _{2}(y_{t}-z_{t}) +\varepsilon _{rt}. \tag{(ref)'} \end{align} Here, only the monetary policy equation is misspecified; other equations remain correct. \\ Issue (i) is illustrated by Table (ref): when the model is misspecified, Kalman filtered shocks can be cross-correlated even for those specified to be independent. In comparison, the Optimal Transport filtered shocks match the model-implied covariance structure. The left panel of Figure (ref) provides further information. Comparing filtered monetary policy shocks $\varepsilon_{rt}$ with the true values, the OTF estimates recover the true series accurately with an $R^2$ of $0.9997$. In contrast, the KF estimates track it less accurately, with an $R^2$ of $0.8035$. \\ \textit{Issue (ii)} is illustrated by the right panel of Figure (ref): the KF fits inflation perfectly with an $R^2$ of $1.0$, whereas the OTF indicates a notable lack of fit with an $R^2$ of $0.5955$. The KF distorts the shocks to achieve this perfect fit. The $R^2$ between true and fitted interest rates are $1.0$ and $0.9966$ for KF and OTF, respectively. The monetary equation is misspecified but the lack of fit mainly appears in the inflation variable. Because this is a simultaneous equation system, misspecification in one equation can affect all variables. In this case the OTF indicates which variables are most affected by the misspecification.\\ Going beyond this simulation illustration, a comparison of latent processes recovered by the KF and OTF using real data for the smets2007 model can be found in Figure (ref) and Supplemental Figure (ref). Moreover, Section (ref) compares actual and OT filtered observables for the LS and SW models and conducts formal specification tests to determine whether the discrepancy between the observed and fitted series is statistically significant. \begin{table}[h] \caption{Misspecified LS model: Covariance Matrices of True and Filtered Shocks} \setlength\tabcolsep{4.0pt} { \begin{tabular}{c|ccc|ccc|ccc} \hline \hline & \multicolumn{3}{c|}{True (unobserved)} & \multicolumn{3}{c|}{Optimal Transport Filter} & \multicolumn{3}{c}{Kalman Filter}\\ \hline & $\varepsilon_r$ & $\varepsilon_g$ & $\varepsilon_z$ & $\varepsilon_r$ & $\varepsilon_g$ & $\varepsilon_z$ & $\varepsilon_r$ & $\varepsilon_g$ & $\varepsilon_z$\\ \hline $\varepsilon_r$ & 0.0802 & 0 & 0 & 0.0802 & 0 & 0 & 0.0582 & 0.0149 & 0.0017\\ $\varepsilon_g$ & 0 & 0.0296 & 0.0674 & 0 & 0.0296 & 0.0674 & 0.0149 & 0.0227 & 0.0688\\ $\varepsilon_z$ & 0 & 0.0674 & 0.5398 & 0 & 0.0674 & 0.5397 & 0.0017 & 0.0688 & 0.5413\\ \hline\hline \end{tabular} } \\ \linespread{0.2} \captionsetup{justification=justified} \caption*{ \textbf{Legend:} True: model implied covariance matrix. Optimal Transport and Kalman Filters: sample covariance matrix of filtered shocks. Simulated sample size $n=10^6$.} \end{table} \begin{figure}[h] \caption{Misspecified LS model: Filtered Shocks and Fitted Variables} \\ \linespread{0.2} \captionsetup{justification=justified} \caption*{ \textbf{Legend:} First 50 observations of the simulated series from model ((ref)) and fitted series using model ((ref)'). Black solid line: True value, Purple dashed line: Optimal Transport Filter, Blue dashed line: Kalman Filter. } \end{figure} \subsection{Discussion: Standard Filtering under Misspecification} Consider a general state-space model as in ((ref)), associated with a state transition density \(p(z_t|z_{t-1})\) and measurement density \(p(y_t|z_{t})\) (omitting $\theta$ for brevity). Standard filtering such as KF and particle filter proceeds recursively as follows: starting at time $t=0$ with initial distribution $p(z_0)$: (a) form beliefs \( p(z_1)=\int p(z_1|z_0)p(z_0)\,dz_0, \) (b) observe $y_1$ and update beliefs using Bayes' rule to obtain the filtered distribution \( p(z_1|y_1)\propto p(y_1|z_1)p(z_1), \) (c) predict \( p(z_2|y_1)=\int p(z_2|z_1)p(z_1|y_1)dz_1, \) (d) continue for $t=2,3,\dots$ However, under misspecification, in step (b) the researcher observes $\tilde{y}_1$ instead of $y_1$, with a predictive density that differs from \( p(y_1|z_1). \) The resulting filtered distribution for $z_1$ is biased because $p({y}_1|z_1)$ is replaced by $p(\tilde{y}_1|z_1)$. The prediction in step (c) is also distorted because it uses this filtered distribution as input. The mismatch occurs at each step when a new observation enters the filter. As a result, the inferred latent variables become difficult to interpret within the structural model. This issue does not improve with a larger samples. To address this mismatch, we propose constructing a model-consistent sample $y_1,\dots,y_n$ using the observed data $\tilde{y}_1,\dots,\tilde{y}_n$. The main idea is that, at each $t$, an optimal transport step is used to map the observation $\tilde{y}_t \sim \tilde{p}(\tilde{y}_t|\tilde{y}_{t-1},\dots)$ into a model-consistent value $y_t \sim p(y_t|y_{t-1},\dots)$ which is closest to the data. Filtering in steps (b) and (c) then proceeds using $y_t$ rather than $\tilde{y}_t$, eliminating the model-data mismatch. The filtered state variables $z_t$ satisfy the model restrictions by construction. Next we focus on linear state-space models to provide a complete analysis of filtering, estimation and inference, forecasting, and model specification testing. Section (ref) then presents a filtering algorithm for nonlinear state-space models based on the same idea outlined above. \begin{comment}In this setting, the data-generating process for $\tilde{y}_t$ is given by the conditional distribution $\tilde{p}(\tilde{y}_t|\tilde{y}_{t-1},\tilde{y}_{t-2},\dots) := \tilde{p}_{t|t-1}(\tilde{y}_t)$, which needs to be estimated. Filtering methods can be used to evaluate the model's predictive distribution $p(y_t|y_{t-1}) := p_{t|t-1}(y_t)$, which is determined by the model and therefore not estimated. \end{comment} \section{{Optimal Transport Filtering, Estimation, and Forecasting}} For the remainder of the paper, $\tilde{y}_t$ denotes the observations collected by the researcher and $y_t$ denotes data generated from ((ref)). The data generating process for $\tilde{y}_t$ is given by the conditional distribution $\tilde{p}(\tilde{y}_t|\tilde{y}_{t-1},\tilde{y}_{t-2},\dots)$, which needs to be estimated. For a parameter value $\theta$, filtering methods can be used to evaluate the model's predictive distribution $p(y_t|y_{t-1},\dots;\theta)$, this is not estimated. To reduce notation, the same $p$ refers to the joint and marginal distributions of $y$ and $z$. \subsection{Optimal Transport Filtering} Consider linear state-space models of the form: \begin{align} y_t = \mu(\theta) + A(\theta) z_t + B(\theta) v_t, \quad z_{t} = C(\theta) z_{t-1} + D(\theta) v_t, \end{align} where $v_t \sim (0,I)$ are white noise, with its higher order moments left unspecified. The dependence on the parameters $\theta$ will be omitted in the Algorithm below to simplify notation. Specification ((ref)) sets a particular linear structure in ((ref)). The number of structural shocks $v_t$ can be greater, equal, or less than the number of observed outcomes $y_t$. Model ((ref)) includes linearized DSGE models and affine term structure models as special cases. The model is characterized by the conditional mean and covariance of $y_t$ and its prediction error at each $t$. To construct the transport map, the proposed filter requires computing their data-implied counterparts, which requires a flexible approximation to the DGP. Under mild conditions (given below), $\tilde{y}_t$ admits an infinite-order vector autoregressive (VAR) representation. A natural auxiliary model is a VAR(k), a sieve approximation of the VAR($\infty$): \begin{align} \tilde{y}_t = \tilde{\mu} + \sum_{j=1}^k \Psi_j [\tilde{y}_{t-j} - \tilde{\mu}] + e_t. \end{align} This is a reduced-form VAR, and $e_t \sim (0,\tilde{\Sigma}_k)$ is a prediction error that may differ in dimension from the structural shocks. As $k \to \infty$, $\tilde{\Sigma}_k$ converges to the innovation covariance matrix of the VAR($\infty$). In practice, the number of lags $k$ should be sufficiently large so that no significant residual autocorrelation remains. See kuersteiner2005, and references therein, for automated lag-length selection procedures. This is a semi-parametric model because $e_t$ is specified up to the first two unconditional moments and is distributionally unrestricted. Algorithm (ref) below presents the OTF, using as inputs: 1) the residuals $\hat{e}_t$ from a vector autoregression of $\tilde{y}_t$ on its $k$ lags, the covariance matrix of these residuals, and 2) for a given parameter value $\theta$, the model-implied state-space coefficients $\mu,A,B,C,D$. Algorithm (ref) only involves matrix operations and can be readily applied to models where the KF is used. \begin{algorithm}[H] \caption{Optimal Transport Filter: Linear State-Space Models} { \begin{algorithmic}[1] \Procedure{\textsc{otf}}\newline \textbf{Inputs:} 1) Sample: data $\tilde{y}_1,\dots,\tilde{y}_n$, residuals $\hat{e}_1,\dots,\hat{e}_n$ from a VAR($k$) regression with $k \geq 1$\newline\hphantom{\textbf{Inputs:}} 2) Model: coefficients $\mu,A,B,C,D$. Initial beliefs $z_0 \sim (m_{0|0},V)$ \newline \textbf{Outputs:} 1) Mapped data $y_1,\dots,y_n$, 2) Filtered states $z_{t|t} \sim (m_{t|t},V)$ \For{$t \in \{1,\dots,n\}$} \State{\textbf{Predict:} $m_{t|t-1} = C m_{t-1|t-1}$, $\mu_{t|t-1} = \mu + A m_{t|t-1}$} \Comment{(KF)} \State{\textbf{Transport:} $y_t = \mu_{t|t-1} + P \hat{e}_t$} \Comment{(OT)}\newline \hphantom{\textbf{Compute:}} where $P = \tilde{\Sigma}_{nk}^{-1/2}[\tilde{\Sigma}_{nk}^{1/2} \Sigma \tilde{\Sigma}_{nk}^{1/2}]^{1/2}\tilde{\Sigma}_{nk}^{-1/2}$ \,\, (Transport Map)\newline \hphantom{\textbf{Compute:}} and $\Sigma = \text{var}_{t|t-1}(y_t)$, $\tilde{\Sigma}_{nk} = \widehat{\text{var}}_{t|t-1}(\tilde{y}_t)$ (Innovation Variance) \State{\textbf{Update:} $m_{t|t} = m_{t|t-1} + K P \hat{e}_t$} \Comment{(KF)}\newline \hphantom{\textbf{Compute:}} where $K = \overline{V} A^\prime \Sigma^\dagger$ (Kalman gain)\newline \hphantom{\textbf{Compute:}} and $\overline{V} = \text{var}_{t|t-1}(z_t)$ \EndFor \EndProcedure \end{algorithmic}} \end{algorithm} Algorithm (ref) combines time-invariant KF iterations with an optimal transport (OT) map $P$. It adjusts the innovations $\hat{e}_t$, whose sample variance is $\tilde{\Sigma}_{nk} = \hat{\text{var}}(\tilde{y}_t|\tilde{y}_{t-1},\dots)$, to match the variance $\Sigma(\theta) = \text{var}(y_t|y_{t-1},\dots)$ implied by model ((ref)). The predictive distributions are summarized by $m_{t|t} = \mathbb{E}(z_t|y_t,\dots,y_1)$, $m_{t|t-1} = \mathbb{E}(z_t|y_{t-1},\dots,y_1)$, $\mu_{t|t-1} = \mathbb{E}(y_t|y_{t-1},\dots,y_1)$, and $V = \text{var}(z_t|y_t,\dots)$; $\tilde{\Sigma}_{nk}^{1/2}$ and $\tilde{\Sigma}_{nk}^{-1/2}$ are the matrix square root of $\tilde{\Sigma}_{nk}$ and its inverse. There are three main steps in the Algorithm: \textbf{Filtering}, \textbf{Transport}, and \textbf{Update}. The following discusses each step in more detail. After that, a simple example will illustrate how the series $y_t$ is constructed and how estimation is performed. \textbf{Predict.} The prediction step is a standard KF operation. The matrices $V,\Sigma$, and $K$ solve the system of equations anderson1979: $\bar{V} = C V C^\prime + DD^\prime,K = \bar{V} A^\prime \Sigma^{\dagger},V = (I - KA)\bar{V},\Sigma = A C V C^\prime A^\prime + (B + AD)(B + AD)^\prime$, where $\Sigma^{\dagger}$ denotes the Moore-Penrose inverse of $\Sigma$ if it is singular, and otherwise its inverse. Besides $V$ and $\Sigma$ defined above, the matrix $\overline{V} = \text{var}(z_t|y_{t-1},\dots)$ measures the one-step-ahead prediction error for $z_t$. These matrices are standard KF quantities. \textbf{Transport.} The transport step maps an observation $\tilde{y}_t$ with conditional distribution $\tilde{p}_{t|t-1}$ to a $y_t$ with conditional distribution $p_{t|t-1}$. This is done by solving the minimization problem: \[ \min_{\pi_{t|t-1}} \mathbb{E}_{\pi_{t|t-1}}(\|\tilde{y}_t - y_t\|^2), \] where ${\pi_{t|t-1}}$ represents any joint distribution with marginal distributions $\tilde{p}_{t|t-1}$ and $p_{t|t-1}$: \[ y_t|\{y_{t-1},y_{t-2},\dots\} \sim ( \mu_{t|t-1}, \Sigma ), \quad \tilde{y}_t|\{\tilde{y}_{t-1},\tilde{y}_{t-2},\dots\} \sim ( \tilde{\mu}_{t|t-1}, \tilde{\Sigma} ),\] where $(\tilde{\mu}_{t|t-1},\mu_{t|t-1})$ are the conditional means of $(\tilde{y}_t,y_t)$ and $(\tilde{\Sigma},\Sigma)$ the forecast error covariance matrices; $(\mu_{t|t-1},\Sigma)$ are produced by the Kalman recursions in the Predict step; and $(\tilde{\mu}_{t|t-1},\tilde{\Sigma})$ are estimated from the data using the vector autoregression. This setup ensures that the generated $y_t$ represents a draw from the model distribution $p_{t|t-1}$. The solution is unique and in closed form, given by $y_t=\mu_{t|t-1}+P\hat{e}_t$, as shown in Step 4 of the Algorithm. All subsequent belief updating and filtering operations are based on the new $y_t$ rather than the original $\tilde{y}_t$, ensuring model consistency. In the absence of misspecification, $P$ is the identity matrix and $y_t$ equals $\tilde{y}_t$. The following provides details on how the transport map is derived and why it is semiparametrically valid. Here, the marginal distributions are only specified up to second moments. As a result, $\pi_{t|t-1}$ are also defined up to second moments: \[ \left( \begin{array}{c} y_t \\ \tilde{y}_t \end{array} \right) \Bigg| \left( \begin{array}{c} y_{t-1},y_{t-2},\dots \\ \tilde{y}_{t-1},\tilde{y}_{t-2},\dots \end{array} \right) \sim \left( \left( \begin{array}{c} \mu_{t|t-1} \\ \tilde{\mu}_{t|t-1} \end{array} \right), \left( \begin{array}{cc} \Sigma & C_{t|t-1} \\ C_{t|t-1}^\prime & \tilde{\Sigma} \end{array} \right) \right), \] where $C_{t|t-1}$ is the conditional covariance between $y_t$ and $\tilde{y}_t$, and $\tilde{\mu}_{t|t-1},{\mu}_{t|t-1},\tilde{\Sigma},$ and $\Sigma$ are the same as above. Recall that the Optimal Transport problem is to minimize $\mathbb{E}_{\pi_{t|t-1}}(\|\tilde{y}_t-y_t\|^2) = \|\tilde{\mu}_{t|t-1} - \mu_{t|t-1}\|^2 + \text{trace}(\tilde{\Sigma}+\Sigma) - 2 \text{trace}(C_{t|t-1})$ in this setup, with the additional constraint that $\pi_{t|t-1}$ is a proper distribution, i.e. $C_{t|t-1}$ cannot be arbitrary. This implies that the optimal transportation problem can be written as a semidefinite program: \begin{align} \min_{C} \Big( -2 \text{trace} ( C ) \Big) \text{ subject to } \left( \begin{array}{cc} \Sigma & C \\ C^\prime & \tilde{\Sigma} \end{array} \right) \geq 0. \end{align}Any transport map between $\tilde{y}_t$ and $y_t$ with covariance $C$ that solves $(\ref{eq:OT_L})$ is optimal. Since the distributions are not fully specified, the map is not uniquely defined. In particular, the linear map which solves the Gaussian case: \begin{align*} T : \tilde{y}_t \to \mu_{t|t-1} + P (\tilde{y}_t - \tilde{\mu}_{t|t-1}), \text{ where } P = \tilde{\Sigma}^{-1/2} \left( \tilde{\Sigma}^{1/2} \Sigma \tilde{\Sigma}^{1/2}\right)^{1/2} \tilde{\Sigma}^{-1/2}, \end{align*} is optimal and preserves the linearity of the process. These derivations follow from dowson1982, olkin1982, and givens1984. See peyre2019, Remark 2.31, for additional discussion of the Gaussian case. The Transport step in Algorithm (ref) uses the plugin estimate $\tilde{\Sigma}_{nk}$ for $P$, and $\hat{e}_t$ for $\tilde{y}_t-\tilde{\mu}_{t|t-1}$. \textbf{Update.} A key difference with KF is in the update step. The standard KF update is $m_{t|t} = m_{t|t-1} + K \tilde{e}_t$ where $\tilde{e}_t$ are prediction errors for $\tilde{y}_t$ computed using model ((ref)) and the $K$ matrix is the Kalman gain. Here, the prediction errors $\hat{e}_t = \tilde{y}_t - \tilde{\mu}_{t|t-1}$ are based on the auxiliary VAR model and are transported using the matrix $P$ to have variance $\Sigma(\theta)$. Enforcing the model-based covariance structure ensures the new data is model-consistent. If the model is misspecified such that $\tilde{e}_t$ are serially correlated and/or do not have covariance matrix $\Sigma$, then the Kalman filtered variable will not satisfy model constraints. \textbf{Accommodating some non-linearities:} The plugin transport map extends to non-linear models of the form: \[ y_t = \mu(x_t;\theta) + \Sigma^{1/2}(x_t;\theta) v_t, \] where $v_t \sim (0,I)$ and $x_t$ is observed, or can be perfectly inferred, at time $t-1$. The solution to ((ref)) now changes with $t$: $P_{t|t-1} = \tilde{\Sigma}_{t|t-1}^{-1/2}( \tilde{\Sigma}_{t|t-1}^{1/2} \Sigma(x_t;\theta) \tilde{\Sigma}_{t|t-1}^{1/2} )^{1/2} \tilde{\Sigma}_{t|t-1}^{-1/2}$; the map becomes $T_{t|t-1} : \tilde{y}_t \to \mu(x_t;\theta) + P_{t|t-1} (\tilde{y}_t - \tilde{\mu}_{t|t-1})$. The requirement that $x_t$ is observable accommodates (G)ARCH but not stochastic volatility models, for instance. Choices of auxiliary models $\tilde{p}_{t|t-1}$ used for these models in simulation-based estimation are referenced below. \subsection{Optimal Transport Estimation (OTE)} An output of the OTF Algorithm (ref) is the model-consistent series $y_t$, which will be referred to as \textit{coupled series}, or \textit{coupling}. The following considers estimating the parameters $\theta$ by minimizing the discrepancy between the original sample $\tilde{y}_t$ and its coupling $y_t$. The coupling $y_t$ depends on two sets of parameters: the structural coefficients $\theta$, and reduced-form auxiliary parameters $\psi_k$. For stationary linear processes approximated by a finite-order VAR(k), $\psi_k = (\tilde{\mu}^\prime,\text{vech}(\tilde{\Sigma})^\prime,\text{vec}(\Psi_1)^\prime,\dots,\text{vec}(\Psi_k)^\prime)^\prime$, where $\text{vec}$, $\text{vech}$ denote the vectorization and half vectorization magnus2019. In practice, OTF relies on OLS estimates $\hat{\psi}_{nk} = (\tilde{\mu}_{n}^\prime,\text{vech}(\tilde{\Sigma}_{nk})^\prime,\text{vec}(\hat{\Psi}_1)^\prime,\dots,\text{vec}(\hat{\Psi}_k)^\prime)^\prime$. The estimation is conducted as follows: given parameters $(\theta,\hat{\psi}_{nk})$, use Algorithm (ref) to generate a series $y_t(\theta;\hat{\psi}_{nk})$ and compute the loss function: \begin{align*} Q_n(\theta;\hat{\psi}_{nk}) = \frac{1}{n} \sum_{t=1}^n \|y_t(\theta;\hat{\psi}_{nk}) - \tilde{y}_t\|^2_{W_n}, \end{align*} for some symmetric positive definite weighting matrix $W_n$. The optimal transport estimator (OTE) is the minimizer $\hat{\theta}_n$ of $Q_n$. For a d-dimensional vector, i.e., $y_t = (y_{t,1},\dots,y_{t,d})$, setting $W_n = \text{diag}(\text{var}(\tilde{y}_{t,1}),\dots,\text{var}(\tilde{y}_{t,d}))^{-1}$ gives the qualitative interpretation that $\hat{\theta}_n$ maximizes the average R-squared between $\tilde{y}_t$ and its coupling $y_t$, i.e. $R^2_j = 1- [\sum_t (y_{t,j} - \tilde{y}_{t,j})^2 ]/[\sum_t (\tilde{y}_{n,j} - \tilde{y}_{t,j})^2 ]$ for $j \in \{1,\dots,d\}$. This choice of $W_n$ was used in all simulated and empirical examples. For DSGE models, it is common to incorporate prior information. This can be accommodated here by penalization: $Q_n(\theta;\hat{\psi}_{nk}) - \frac{1}{n} \log(\pi(\theta))$, where $\pi$ is the prior density. Under suitable regularity conditions, the first-order asymptotic properties of $\hat{\theta}_n$ are unchanged. Although the OTE might appear abstract, Theorem (ref) in Section (ref) implies that under regularity conditions, given below, the sample loss $Q_n$ converges in probability to: \begin{align*} Q(\theta;\psi_0) = \|\tilde{\mu} - \mu(\theta)\|^2_W + \sum_{j=0}^\infty \text{trace} \left( \tilde{\Sigma}^{1/2} \{ \tilde{\Lambda}_j - \Lambda_j(\theta) P(\theta;\tilde{\Sigma}) \}^\prime W \{ \tilde{\Lambda}_j - \Lambda_j(\theta)P(\theta;\tilde{\Sigma}) \} \tilde{\Sigma}^{1/2} \right), \end{align*} where $W_n$ converges to $W$; $\tilde{\Lambda}_0=\Lambda(\theta)_0=I_d$; $(\tilde{\Lambda}_j,\Lambda(\theta)_j)_{j\geq 1}$ are the vector moving average coefficients; $(\tilde{\mu},\mu(\theta))$ the unconditional means; $(\tilde{\Sigma},\Sigma(\theta))$ the forecast error covariances of $(\tilde{y}_t,y_t(\theta))$; and $P(\theta;\tilde{\Sigma})$ is the transportation matrix introduced above. From this expression, the OTE can be interpreted as minimizing a distance between the sample and model's impulse responses at all horizons. This is illustrated with a simple example below. The OTE is closely related to simulation-based estimators: the simulated method of moments, indirect inference, and the efficient method of moments. Both rely on the original sample $\tilde{y}_t$ and model-implied data $y_t(\theta)$. An important difference is that OT constructs $y_t$ from reduced-form innovations $\hat{e}_t$ evaluated from the sample $\tilde{y}_t$, whereas simulation-based estimator draw structural innovations from a known parametric distribution. Similar to these methods, OT estimation could minimize a distance between moments computed from $\tilde{y}_t$ and $y_t(\theta)$, respectively. A common criticism about fitting a particular set of moments is that the corresponding coefficients may not provide a good fit on non-fitted moments. The following Lemma addresses this particular concern. \begin{lemma} Suppose $\tilde{y}_t$ and $y_t(\theta;\psi_0)$ are covariance stationary. Let ${\bf{\tilde{y}}}_t = (\tilde{y}_t,\tilde{y}_{t-1},\dots)$, ${\bf{y}}_t(\theta;\psi_0) = (y_t(\theta;\psi_0),y_{t-1}(\theta;\psi_0),\dots)$ denote the infinite history of $\tilde{y}_t$ and ${y}_t(\theta;\psi_0)$. Let: \begin{align*} \mathcal{M} = \Big\{ m:{\bf{\tilde{y}}}_t \to \sum_{j \geq 0} A_{j}^\prime \tilde{y}_{t-j} + \sum_{j,j^\prime \geq 0} \tilde{y}_{t-j}^\prime B_{j,j^\prime} \tilde{y}_{t-j^\prime},\text{ s.t. } \sum_{j=0}^\infty \|A_j\|^2 \leq M_1^2, \sum_{j,j^\prime \geq 0} \|B_{j,j^\prime}\|_{\infty} \leq M_2 \Big\}, \end{align*} denote the set of quadratic moments over the infinite histories. Then: \[ \sup_{m \in \mathcal{M}}\Big\|\mathbb{E}\left( m({\bf{\tilde{y}}}_t)\right) - \mathbb{E} \left( m({\bf{{y}}}_t(\theta;\psi_0)) \right) \Big\| \leq \sqrt{C_W} (M_1+M_2 \sqrt{\mathbb{E}(\|\tilde{y}_t\|^2)}) \sqrt{Q(\theta;\psi_0) } + C_W M_2 Q(\theta;\psi_0), \] for some $C_W$ which depends on $W$. If $\tilde{\mu} = \mu(\theta)$, $M_1$ can be removed from the upper bound. \end{lemma} Lemma (ref) shows that the OT loss, $Q$, bounds the worst-case fit over a large class of quadratic moments which includes means, covariances, and autocovariances. These are routinely used in (simulated) minimum distance (SMD) estimations of ((ref)), where the fit is not guaranteed for moments that do not enter the estimation. Similarly, the Kullback-Leibler divergence, used in MLE, does not provide this type of guarantee for moments. \textbf{Accommodating some non-linearities (cont'd):} For non-linear models discussed in the previous subsection, the choice of $\tilde{p}$ is related to the choice of auxiliary models in simulation-based estimation, and must provide an adequate approximation to the DGP. gallant1996 suggest several models, including the SNP estimator of gallant1987 and a nonparametric ARCH model.\footnote{gallant1987 call $\tilde{p}$ the score generator for the Efficient Method of Moments.} altissimo2009, kristensen2012 consider kernel-density estimates when the model is Markovian in the observables. Once chosen, their parameter estimates represent $\hat{\psi}_{nk}$, and the estimation of structural parameters proceeds in the same way, using the same objective function as in the linear case. \subsection{A Pen & Pencil Example} The following example illustrates how the OTF constructs a model-consistent series, how the parameters are estimated, and how the estimates compare with quasi-maximum likelihood. Consider a DGP and a structural model given, respectively, by \begin{align*} \tilde{y}_t &= \tilde{\mu} + \rho \tilde{y}_{t-1} + \tilde{\sigma} e_t, \qquad y_t = \mu + \sigma \left( v_t + \beta_1 v_{t-1} + \dots + \beta_q v_{t-q} \right), \end{align*} where $e_t \sim (0,1)$ and $v_t \sim (0,1)$. An MA($q$) model is thus fitted to data generated by an AR(1). Assume the lag order $q$ is finite, so that the model is misspecified. Also assume the MA($q$) parameters specify a stationary and invertible process so that $v_t$ are the prediction errors, that is, the shocks in the Wold decomposition. Use an AR($k$) model with $k \ge 1$ as the auxiliary model, and let $\hat{e}_t$ denote the standardized residuals from this regression. By Algorithm (ref), for a given value of $\theta = (\mu,\beta_1,\dots,\beta_q,\sigma)$, the OTF maps the prediction errors from the AR($k$) into the prediction errors of the MA($q$). This yields the following coupled series: \( y_t = \mu + \sigma \left[\hat{e}_t + \beta_1 \hat{e}_{t-1} + \dots + \beta_q \hat{e}_{t-q}\right]. \) Accordingly, the OTE for $\theta$ minimizes the loss function \[ Q_n(\theta;\hat{\psi}_{nk}) = \frac{1}{n}\sum_{t=1}^n (\tilde{y}_t-y_t)^2 = \frac{1}{n}\sum_{t=1}^n (\tilde{y}_t - \mu - \sigma \hat{e}_t - \sigma \beta_1 \hat{e}_{t-1} - \dots - \sigma \beta_q \hat{e}_{t-q})^2, \] which amounts to regressing $\tilde{y}_t$ on $(\hat{e}_t,\hat{e}_{t-1},\dots,\hat{e}_{t-q})$ and an intercept. Since the difference between $(\hat{e}_t,\dots,\hat{e}_{t-q})$ and $(e_t,\dots,e_{t-q})$ is asymptotically negligible, $\hat{\theta}_n$ converges in probability to $\theta_0 = (\tilde{\mu},\rho,\rho^2,\dots,\rho^q,\tilde{\sigma})$. The OTE recovers the impulse response function of $\tilde{y}_t$ to $e_t$ up to the specified order $q$, even though the estimated model is misspecified. The maximum likelihood estimator of $\theta$ does not have this property. For example, when $q=1$, conditional on $v_{0} = 0$ the maximum likelihood estimator minimizes: \[ L_n(\theta) = \frac{1}{n} \sum_{t=1}^n ( \tilde{y}_t - \mu + \sum_{k=1}^{t-1} (-\beta_1)^k \tilde{y}_{t-k} )^2/\sigma^2 + n \log(\sigma^2), \] which does not have closed-form expressions for $\beta_1$, $\sigma^2$ when $\rho \neq 0$. \begin{figure}[h] \caption{Misspecified MA(1) model: MLE and OTE estimates ($n=10^4$)} \\ \linespread{0.2} \captionsetup{justification=justified} \caption*{ \textbf{Legend:} Estimates from a simulated sample of $n=10^4$ observations, with $q=1$. Black solid line: MLE estimates, Purple solid line: Optimal Transport Estimates (OTE). } \end{figure} Figure (ref) compares MLE and OTE for an MA(1) model as a function of the true AR(1) coefficient $\rho \in [0,0.95]$, fixing $(\tilde{\mu},\tilde{\sigma})=(0,1)$. As expected, the OTE of the MA(1) coefficient $\beta_1$ varies linearly with $\rho$, since it converges to $\rho$ asymptotically. It recovers the impulse response at horizon $h=1$. The MLE of $\beta_1$ depends nonlinearly on $\rho$ and systematically underestimates the impulse response at horizon $h=1$. For $\sigma$, the OT estimate matches the true value $\tilde{\sigma}=1$, whereas the ML estimate diverges nonlinearly as $\rho$ increases. \subsection{Optimal Transport Forecasting} The standard practice to forecast from a sample ($\tilde{y}_1,\dots,\tilde{y}_n$) using model ((ref)) is to apply a standard filter to recover latent variables $(z_1,\dots,z_n)$, and use the model dynamics to predict $z_{n+h|n}(\theta)$ and then $y_{n+h|n}(\theta)$ at $h \geq 1$. As shown above, when the true data generating process does not coincide with ((ref)), the filtered variables are not model-consistent and, by extension, neither are the resulting forecasts. Under misspecification, mismatch between model and data implies that the resulting forecast cannot be interpreted through the structural model. Using Algorithm 1, there are \textit{a priori} two approaches to constructing a model-consistent forecast conditional on the parameter estimates: \textbf{(i) direct forecast}: predict $y_{n+h|n}^{D}(\theta;\hat{\psi}_{nk})$ using the model-implied conditional mean of the observable in model ((ref)), conditional on the generated sample $y_1(\theta;\hat{\psi}_{nk}),\dots,y_n(\theta;\hat{\psi}_{nk})$; \textbf{(ii) indirect forecast}: predict $\tilde{y}_{n+h|n}$ using the conditional mean given by the auxiliary model evaluated at $\hat{\psi}_{nk}$, and then apply Algorithm 1 to transport the extended sample $\tilde{y}_1,\dots,\tilde{y}_n,\tilde{y}_{n+1|n},\tilde{y}_{n+2|n},\dots$ to produce the forecast $y_{n+h|n}^{D}(\theta;\hat{\psi}_{nk})$. The following lemma shows that the two forecasts are identical. \begin{lemma} For any $\theta \in \Theta$ and $h\geq 1$, the following holds: $y_{n+h|n}^{I}(\theta;\hat{\psi}_{nk}) = y_{n+h|n}^{D}(\theta;\hat{\psi}_{nk}).$ \end{lemma} An implication of Lemma (ref) is that the reduced form forecast can be decomposed into $\tilde{y}_{n+h|n} = y_{n+h|n}^{D} + u_{n+h|n}$ where $y_{n+h|n}^{D}$ decomposes the forecast in terms of strutural quantities and $u_{n+h|n}$ is the variation that the structural model cannot explain. Among competiting structural models, those with a larger $R^2$ can explain a greater share of the reduced-form forecast. \section{Related Literatures} Textbook references on optimal transport (OT) include villani2003 for theory, peyre2019 for computation, and galichon2018 for Economics. In statistics, much of the methodology and theory considers OT between iid samples. dudley1969 showed that the empirical Wasserstein distance suffers from a curse of dimensionality, unlike the plug-in approach used here. The literature is much more limited for dependent data. o2022 construct couplings between finite state Markov Chains using dynamic programming methods. They do not consider parameter estimation. Several papers consider parameter estimation using the Wasserstein distance. bassetti2006a and bassetti2006b study the estimation of location and scale for univariate distributions. They derive consistency and a non-standard limiting distribution for the estimator; bernton2019a extend their results. As a minimum-distance estimator, alternatives to OTE include the Simulated Method of Moments, Indirect Inference gourieroux1996, and adversarial estimation using GANs kaji2023. genevay2018 discuss the advantages of using OT over classifiers found in GANs. forneron2023 considers semi-nonparametric simulation-based estimation, but assumes correctly specified dynamics. These methods do not recover the latent variables. Algorithm (ref) is closely related to a goodness-of-fit plot in the Real Business Cycle literature. plosser1989 and king1999 compute historical productivity shocks outside the model and use them to simulate a one-shock RBC economy. They plot simulated against real data to show the fit of calibrated models. Misspecification robust filtering also considers model misspecification but aims to recover the true latent variable. The main goal is to reduce sensitivity to local misspecification over a pre-specified neighborhood, see e.g. sayed2001 and shafieezadeh2018. This relates to hansen2008's approach to robustness in Economics. Here, the model can be globally misspecified; the filtered values are computed under model constraints. There is a rich literature on quasi-maximum likelihood estimation under misspecification starting with white1982. Result on quasi-ML estimation of hidden Markov models are fewer: mevel2004 consider discrete state-space models, douc2012 give conditions for parameter consistency for more general state-space models. These estimators rely on Kalman or particle filters which, as illustrated above, recover latent variables that typicall do not satisfy moment contraints under misspecification. There a also a number of misspecification-robust moment-based estimators, e.g. schennach2007 or antoine2021 which fit to model to a selected set of moments. Several papers consider estimation and policy analysis with misspecified DGSE models. del2007 and del2009 use a DSGE-VAR framework where a hyperparameter penalizes between a reduced form and structural model. Here, the flexible VAR is used to enforce the model structure with the coupling. This ensures the parameters are internally valid, i.e. characterize the dynamics of $y_t$. watson1993 proposed a goodness of fit measure based on the magnitude (in the frequency domain) of the measurement error required for the model to match the second-order moments of the data. He does not consider filtering, estimation, inference or specification testing. \section{Large Sample Properties of OTE} The following derives consistency and asymptotic normality results for a class of linear processes, which includes linear state-space models described by ((ref)): \begin{align} y_t = \mu(\theta) + A(\theta) z_t + B(\theta) v_t, \quad z_{t} = C(\theta) z_{t-1} + D(\theta) v_t. \tag{(ref)} \end{align} \paragraph{Notation:} The parameters $\theta \in \Theta \subseteq \mathbb{R}^{d_\theta}$. Norms: for a matrix $A = (a_{ij})$ of size $n \times m$, the baseline norm is $\|A\| = \sqrt{\text{trace}(A'A)}$, the operator norm is $\|A\|_{op} = \sqrt{\lambda_{\max}(A'A)}$, the sup norm is $\|A\|_{\infty} = \max_{i,j} |a_{ij}|$. Eigenvalues: for a symmetric matrix $A$ of size $n \times n$, $\lambda_j(A)$ denotes the $j$-th eigenvalue, $1 \leq j \leq n$, in increasing order, $\lambda_{\max}(A) = \lambda_n(A)$ and $\lambda_{\min}(A) = \lambda_1(A)$; $\underline{\lambda} \preceq A$ implies $\lambda_{\min}(A) \geq \underline{\lambda}$ and $A \preceq \overline{\lambda}$ implies $\lambda_{\max}(A) \leq \overline{\lambda}$. For a matrix $A$ of size $n \times m$, the singular values are given by $\sigma_j(A) = \sqrt{\lambda_{j}(A'A)}$ if $m < n$, or $\sigma_j(A) = \sqrt{\lambda_{j}(AA')}$ if $n < m$, for $j = 1,\dots,\min(n,m)$; $\sigma_{\min}(A) = \sigma_1(A)$ and $\sigma_{\max}(A) = \sigma_{\min(m,n)}(A) = \|A\|_{op}$. \subsection{Consistency and Asymptotic Normality} The true data-generating process (DGP) and the model are assumed to be left-invertible; i.e. they admit a one-sided infinite vector moving-average (VMA) representation.\footnote{Some models can feature non-invertibility, this is the case with permanent income fernandez2007. Algorithm (ref) only relies on second-order moments whereas identification and estimation of non-invertible models rely on higher-order cumulants, which is beyond the scope of this paper.} \begin{assumption} $\tilde{y}_t$ and $y_t$ admit causal VMA$(\infty)$ representations: \begin{align*} \tilde{y}_t = \tilde{\mu} + e_t + \sum_{j=1}^\infty \tilde{\Lambda}_j e_{t-j}, \quad y_t(\theta) = \mu(\theta) + \xi_t + \sum_{j=1}^\infty \Lambda_j(\theta) \xi_{t-j}, \end{align*} for any $\theta \in \Theta$, where $e_t$ and $\xi_t$ are white noise with variance $\tilde{\Sigma}$ and $\Sigma(\theta)$. \end{assumption} The VMA innovations $\xi_t$ need not coincide with, or span, the structural innovations $v_t$ fernandez2007. The number of structural innovations can be greater than, equal to, or less than the number of observables. In the latter case, the models are stochastically singular.\footnote{For instance: multivariate RBC models with a single shock to productivity are stochastically singular. See komunjer2011, qu2018 for identification and estimation with stochastic singularity.} Our empirical applications cover all three situations. Using the VMA representation, Algorithm (ref) involves the following quantities: $\tilde{\Sigma} = \text{var}(e_t)$, $\Sigma(\theta) = \text{var}(\xi_t)$, $\tilde{\mu}_{t|t-1} = \tilde{\mu} + \sum_{j=1}^\infty \tilde{\Lambda}_j e_{t-j}$, and $\mu_{t|t-1} = \mu(\theta) + \sum_{j=1}^\infty \Lambda_j(\theta) \xi_{t-j}$. Note that $\Sigma(\theta)$ is the same as that in Algorithm (ref), where $y_t$ admits a state-space representation. Take the transport map $P(\theta;\tilde{\Sigma})$, computed using $\tilde{\Sigma}$ and $\Sigma(\theta)$, the coupled series $y_t$ is: \[ y_t(\theta;\psi_0) = \mu(\theta) + P(\theta;\tilde{\Sigma}) e_t + \sum_{j=1}^\infty \Lambda_j(\theta) P(\theta;\tilde{\Sigma}) e_{t-j}. \] The index $\psi_0$ refers to the innovations and the variance used to compute the coupling. For $\psi = \psi_0$, the true errors $e_t$ and $\tilde{\Sigma}$ are used, as above. For $\psi = \hat{\psi}_{nk}$, the residuals $\hat{e}_t$ and sample variance $\tilde{\Sigma}_{nk}$ are used, with the convention that $\hat{e}_t = 0$ for $t \leq 0$. For $\psi = \psi_{k}$, the error $e_{t,k}$ and $\tilde{\Sigma}$ are used, where $e_{t,k} = \tilde{y}_t - \tilde{\mu} - \sum_{j=1}^k \tilde{\Psi}_j[\tilde{y}_{t-j}-\tilde{\mu}]$ are the VAR($k$) errors. The KF steps in Algorithm (ref) compute the VMA coefficients $\Lambda_j(\theta)$ using the state-space representation ((ref)). Suppose $\|C(\theta)\|_{op} < 1$, iterate on the KF and OT steps to find: \[ y_t(\theta;\psi_0) = \mu(\theta) + P(\theta;\tilde{\Sigma}) e_t + \sum_{j=1}^\infty A(\theta)C^j(\theta) K(\theta) P(\theta;\tilde{\Sigma}) e_{t-j}, \] where $K(\theta)$ is the Kalman gain. Thus, $\Lambda_0 = I_d$ and $\Lambda_j(\theta) = A(\theta)C^j(\theta) K(\theta)$ for each $j \geq 1$. \begin{assumption} (i). $\sum_{j=1}^\infty j^{1/2}\|\tilde{\Lambda}_j\|<\infty$ and $\text{det}\left( \sum_{j=0}^\infty \tilde{\Lambda}_j z^j \right) \neq 0$ for all $|z| \leq 1$ with $z \in \mathbb{C}$; (ii). $e_t$ is strictly stationary, $\mathbb{E}_{t-1}(e_t)=0$, $\mathbb{E}(e_te_t^\prime) = \tilde{\Sigma}$, and $0 < \underline{\lambda} \preceq \tilde{\Sigma} \preceq \overline{\lambda} < \infty$; (iii). for some $r>4$, $\mathbb{E}(\|e_t\|^{2r}) < \infty$, and $e_t$ is $\alpha$-mixing with size $-a$, where $a > r/(r-2)$. \end{assumption} Assumption (ref) provides several sufficient conditions for $\tilde{y}_t$ to admit a VAR($\infty$) representation and to study the OLS estimates hannan2012. The mixing conditions are needed to derive near-epoch dependence (NED) properties for $y_t(\theta;\hat{\psi}_{nk})$, its derivatives, and asymptotic results for $\tilde{\Sigma}_{nk}$. Assumption (ref) allows for unmodelled dependence in higher-order moments, such as conditional heteroskedasticity (ARCH, GARCH) or stochastic volatility that satisfy a strong-mixing condition. \begin{assumption} $\Theta$ is convex and compact and $\theta \to (\mu(\theta),\Sigma(\theta),\Lambda_1(\theta),\dots)$ is three times continuously differentiable, such that: (i). $\text{rank}[\Sigma(\theta)] = r_\Sigma$ for all $\theta \in \Theta$, and $0 \preceq \Sigma(\theta) \preceq \overline{\lambda} < \infty$; (ii). $\sup_{\theta \in \Theta} \|\mu(\theta)\| < \infty$ and $\sum_{j=0}^\infty \sup_{\theta \in \Theta} \|\Lambda_j(\theta)\|_{op} < \infty$; (iii). for $s=1,\dots,3$ and any $i_1,\dots,i_s \in \{1,\dots,d_\theta\}$, $\sup_{\theta \in \Theta} \|\partial^s_{\theta_{i_1},\dots,\theta_{i_s}} \mu(\theta)\|< \infty$, $\sup_{\theta \in \Theta} \| \partial^s_{\theta_{i_1},\dots,\theta_{i_s}} \text{vec}[\Sigma(\theta)]\|_{\infty}< \infty$, and $ \sum_{j=0}^\infty \sup_{\theta \in \Theta} \| \partial^s_{\theta_{i_1},\dots,\theta_{i_s}} \text{vec}[\Lambda_j(\theta)]\|_{\infty} < \infty$. \end{assumption} \begin{assumption} There exists $C \geq 0$, $b \geq 2$, and $\varepsilon > 0$ such that for $s=1,\dots,3$ and any $i_1,\dots,i_s \in \{1,\dots,d_\theta\}$: $\sum_{j = m + 1}^\infty \sup_{\theta \in \Theta} \|\Lambda_j(\theta)\|_{op} \leq C m^{-(b+\varepsilon)}$, $\sum_{j = m + 1}^\infty \|\tilde{\Lambda}_j\|_{op} \leq C m^{-(b+\varepsilon)}$, and $\sum_{j = m + 1}^\infty \sup_{\theta \in \Theta} \|\partial^s_{\theta_{i_1},\dots,\theta_{i_s}} \text{vec}[\Lambda_j(\theta)]\|_{\infty} \leq C m^{-(b+\varepsilon)}$. \end{assumption} Assumptions (ref) and (ref) restrict the dependence of $\tilde{y}_t$, $y_t$, and its derivatives. The constant rank condition is discussed below. The following Lemma gives conditions on the state-space representation ((ref)) for which Assumption (ref) holds. Lemma (ref) in Appendix (ref) further shows that Assumption (ref) also holds for any $b \geq 2$ and $\varepsilon > 0$, with an appropriate constant $C > 0$. \begin{lemma}[State-Space Model - VMA representation] If $\Theta$ is convex and compact, and the following conditions hold: (i). $\text{rank}[\Sigma(\theta)] = r_{\Sigma}$ for all $\theta \in \Theta$; (ii). $\Sigma(\cdot)$, $\mu(\cdot)$, $A(\cdot)$, $B(\cdot)$, $C(\cdot)$, and $D(\cdot)$ are three times continuously differentiable with bounded derivatives; (iii). $\inf_{\theta \in \Theta} \inf_{|z| \leq 1, z \in \mathbb{C}}|\text{det}(I - C(\theta)z )| > 0$, then Assumption (ref) holds. \end{lemma} The constant rank condition $\text{rank}[\Sigma(\theta)]=r_{\Sigma}$, which appears in Assumption (ref) and Lemma (ref), and the full rank condition $0 < \underline{\lambda} \preceq \tilde{\Sigma}$ (Assumption (ref)) are particularly important for the transport map to be well behaved. Much like the square root of a scalar, $x \to \sqrt{x}$, the matrix square root $A \to A^{1/2}$, used in the transport map, is not continuously differentiable at a singular $A$. The following Lemma derives a new result for the differentiability of $\theta \to A(\theta)^{1/2}$ when $A(\cdot)$ is singular with a constant rank. The proof involves a constructive local block decomposition which can be used to compute differentials analytically. Unlike an eigenvalue decomposition, the block decomposition is smooth under multiplicity of eigenvalues. This should be of independent interest as the matrix square root appears in a variety of settings. \begin{lemma}[Matrix Square Root, Constant Rank] Suppose $\Theta \subset \mathbb{R}^{d_\theta}$ is convex and compact and $\theta \to A(\theta) \geq 0$ is $s$-times continuously differentiable for some $s \geq 1$. Assume that $A(\theta)$ has constant rank $r$, where $1 \leq r \leq d = \text{dim}(A)$ and $0 < \underline{\lambda} \leq \inf_{\theta} \lambda_{r}[A(\theta)] \leq \sup_{\theta} \lambda_{\max}[A(\theta)] \leq \overline{\lambda} < \infty$. Then: (i). There exists $\delta > 0$, such that for any $\theta_0 \in \Theta$, there exists $M(\theta)$ and $B(\theta)$ that are $s$-times continuously differentiable on $\mathcal{B}_\delta(\theta_0) = \{\theta \in \Theta, \|\theta-\theta_0\| \leq \delta\}$, such that $0 < \underline{\lambda}_B \preceq B(\theta) \preceq \overline{\lambda}_B < \infty$, $M(\theta)M(\theta)^\prime = I_d$, and $A(\theta) = M(\theta) \text{blockdiag}[B(\theta), 0_{m,m}] M(\theta)^\prime$ where $m = d-r$. (ii). For all $\theta_0 \in \Theta$, the square root $A(\theta)^{1/2} = M(\theta) \text{blockdiag}[B(\theta)^{1/2}, 0_{m,m}] M(\theta)^\prime$ is $s$-times continuously differentiable on $\mathcal{B}_\delta(\theta_0)$. (iii). The square root $\theta \to A(\theta)^{1/2}$ is $s$-times continuously differentiable on $\Theta$. \end{lemma} The KF recursions are well defined under stochastic singularity anderson1979; however, the likelihood is not defined. If the constant rank condition fails, the transport map becomes non-smooth and the KF steps in Algorithm (ref) become sensitive to numerical accuracy and can be unstable anderson1979. \begin{lemma}[Data: VAR($\infty$) representation, VAR($k$) approximation] Suppose Assumptions (ref), (ref), and (ref) hold. Then $\tilde{y}_t$ admits a VAR($\infty$) representation: $\tilde{y}_t = \tilde{\mu} + \sum_{j=1}^\infty \Psi_j (\tilde{y}_{t-j} - \tilde{\mu}) + e_t$, where $\sum_{j=1}^\infty j^{1/2}\|\Psi_j\|<\infty$, $\Psi_0 = I_d$, and $\text{det}\left( \sum_{j=0}^\infty \Psi_j z^j \right) \neq 0$ for any $|z| \leq 1$. Further, suppose $k \to \infty$ such that $k^3/n \to 0$ and $\sqrt{n} \sum_{j = k+ 1}^\infty \|\Psi_j\| \to 0$, and let $\tilde{\Sigma}_{nk} = \frac{1}{n} \sum_{t=1}^n \hat{e}_t \hat{e}_t^\prime$ and $\tilde{\Sigma}_{n} = \frac{1}{n} \sum_{t=1}^n e_t e_t^\prime$. Then: (i). $\max_{j=1,\dots,k} \|\hat{\Psi}_j - \Psi_j\| = O_p( \sqrt{\log(n)/n})$; (ii). $\tilde{\Sigma}_{nk} - \tilde{\Sigma}_{n} = o_p(1/\sqrt{n})$; (iii). $\tilde{y}_n - \tilde{\mu} = O_p(n^{-1/2})$ and $\tilde{\Sigma}_n - \tilde{\Sigma} = O_p(n^{-1/2})$. \end{lemma} Lemma (ref) combines several existing results for the auxiliary parameters $\hat{\psi}_{nk}$ from the literature, mainly lewis1985 and hannan2012. The conditions on the order of the VAR order, $k$, depend on the decay of the VAR coefficients. If model ((ref)) is correctly specified and the conditions for Lemma (ref) hold - or if the true model is a finite order stationary VARMA - then $\|\Psi_j\|_{op} = O(\bar{\rho}^j)$ for some $\bar{\rho} \in [0,1)$ and $\sqrt{n} \sum_{j=k+1}^\infty \|\Psi_j\| = o(1)$ as long as $\log(n)/k \to 0$. In these cases, the order $k$ can increase very slowly. \begin{theorem}[Consistency] Suppose Assumptions (ref)-(ref) hold, $k$ satisfies the conditions of Lemma (ref), $W_n \overset{p}{\to} W > 0$, and $Q(\theta;\psi_0) = \lim_{n \to \infty} \frac{1}{n} \sum_{t=1}^n \mathbb{E}\left( \|y_t(\theta;\psi_{0}) - \tilde{y}_t \|_{W}^2 \right)$ is uniquely minimized at $\theta = \theta_0$. If $k$ is such that $\sqrt{n}\sum_{j=k+1}^\infty \sup_{\theta \in \Theta} \|\Lambda_j(\theta)\|_{op} = o(1)$, then $\hat{\theta}_n \overset{p}{\to} \theta_0$. \end{theorem} Theorem (ref) shows that the estimator $\hat{\theta}_n$ is consistent for the minimizer $\theta_0$ of $Q(\cdot;\psi_0)$. \begin{theorem}[Asymptotic Normality] Suppose the conditions for Theorem (ref) hold with $\theta_0 \in \text{interior}(\Theta)$. Let $u_{t,k} = y_t(\theta_0;\psi_k) - \tilde{y}_t$, $u_t = y_t(\theta_0;\psi_0) - \tilde{y}_t$, $G_t(\theta_0;\psi_k) = \text{vec}[ \partial_\theta y_t(\theta_0;\psi_k)^\prime ]$ and \begin{align*} M &= \mathbb{E}\left( \partial_\theta y_t(\theta_0;\psi_0)^\prime W \partial_\theta y_t(\theta_0;\psi_0) \right) + \mathbb{E}\left( \left[ u_{t}^\prime W \otimes I \right]\partial_\theta G_t(\theta_0;\psi_0) \right),\\ D_{\theta,\psi}(k) &= \mathbb{E}\left[ \partial_\theta y_t(\theta_0;\psi_k)^\prime W \partial_\psi y_t(\theta_0;\psi_{k0}) + \left[ u_{t,k}^\prime W \otimes I \right]\partial_\psi G_t(\theta_0;\psi_{k0}) \right]. \end{align*} Suppose $M$ is invertible and there exists $\underline{k} \geq 1$ and $c_1 > 0$ such that for all $k \geq \underline{k}$: $0 < c_1 \leq \sigma_{\min}[D_{\theta,\psi}(k)] < \infty$. Define $Z_{k,t} = ((\tilde{y}_t-\tilde{\mu})^\prime,\text{vec}[ e_t \tilde{Y}_{t-1,k}^\prime\Gamma_k^{-1} ]^\prime,\text{vech}[e_te_t^\prime - \tilde{\Sigma}]^\prime)^\prime$, with $\tilde{Y}_{t-1,k} = ((\tilde{y}_{t-1}-\tilde{\mu})^\prime,\dots,(\tilde{y}_{t-k}-\tilde{\mu})^\prime)^\prime$ and $\Gamma_k = \mathbb{E}( \tilde{Y}_{t-1,k}\tilde{Y}_{t-1,k}^\prime)$. Then, the sequence of covariance matrices $V_{n,k} = M^{-1} \text{var} \left[ \frac{1}{\sqrt{n}} \sum_{t=1}^n \{\partial_\theta y_t(\theta_0;\psi_k)^\prime W u_{t,k} + D_{\theta,\psi}(k)Z_{k,t} \}\right] M^{-1},$ is bounded from above. If, in addition, $V_{n,k}^{-1} = O(1)$, then: $\sqrt{n}V_{n,k}^{-1/2}(\hat{\theta}_n - \theta_0) \overset{d}{\to} \mathcal{N}(0,I).$ \end{theorem} Theorem (ref) establishes the asymptotic normality of the estimates $\hat{\theta}_n$. The invertibility of $M$ and the lower bound $c_1$ are local identification conditions. The requirement $V_{n,k}^{-1} = O(1)$ is standard for central limit theorems white2014. The boundedness of $V_{n,k}$ implies a $\sqrt{n}$-rate of convergence for $\hat{\theta}_n - \theta_0$. Note that this rate does not apply to all functionals of $\hat{\psi}_{nk}$; some may converge more slowly lewis1985. The residual $u_t = y_t(\theta_0;\psi_0) - \tilde{y}_t$ measures the model-data discrepancy. Like in OLS, it reflects a simple decomposition of $\tilde{y}_t$ into fitted values $y_t(\theta_0;\psi_0)$ and residuals $u_t$. The $R^2$ introduced earlier measure their relative magnitudes. With $n = \infty$, $R^2 = 1$ indicates correct specification. Formal specification testing is considered in the next subsection. \paragraph{Computing standard errors.} The following describes how to compute standard errors assuming correct specification and allowing for misspecification. For models considered in Section (ref), bootstrap inference could be rather computationally cumbersome. The plugin estimates for $M,D_{\theta,\psi}(k)$, etc, are shown to be consistent in the proof of Theorem (ref). Consistency of standard errors can be deduced from these results, the proof is omitted for brievety. Under correct specification, $u_t = 0$ and $u_{t,k} = o(n^{-1/2})$ don't contribute to the asymptotic standard errors. First, evaluate $\partial_{\theta} y_t(\hat{\theta}_n;\hat{\psi}_{nk})$ and $\partial_{\psi} y_t(\hat{\theta}_n;\hat{\psi}_{nk})$ with finite-differences or by automatic differentiation. Then, compute $\hat{M}_{n} = \frac{1}{n} \sum_{t=1}^n \partial_{\theta} y_t(\hat{\theta}_n;\hat{\psi}_{nk})^\prime W_n \partial_{\theta} y_t(\hat{\theta}_n;\hat{\psi}_{nk})$ and $\hat{D}_{n,\theta,\psi}(k) = \frac{1}{n} \sum_{t=1}^n \partial_{\theta} y_t(\hat{\theta}_n;\hat{\psi}_{nk})^\prime W_n \partial_{\psi} y_t(\hat{\theta}_n;\hat{\psi}_{nk})$. Let $L_n(\hat{\psi}_{nk})$ denote the Gaussian quasi-likelihood for the auxiliary VAR($k$) model, $H_n(\hat{\psi}_{nk})$ its Hessian, and $\partial_{\psi} L_t(\hat{\psi}_{nk})$ the score for $\tilde{y}_t$. Take $\hat{S}_{t,k} = \hat{M}_{n}^{-1} \hat{D}_{n,\theta,\psi}(k) H_n(\hat{\psi}_{nk})^{-1} \partial_{\psi} L_t(\hat{\psi}_{nk})$. Then, $\hat{V}_{n,k}$ is the HAC estimator for the long-run variance of $\hat{S}_{t,k}$. Standard errors are computed from $\hat{V}_{n,k}/n$ in a standard fashion. Allowing for misspecification requires estimating several additional terms. Compute $\hat{u}_t = \hat{u}_{t,k} = y_t(\hat{\theta}_n;\hat{\psi}_{nk}) - \tilde{y}_t$, $\partial_\theta \hat{G}_t(\hat{\theta}_n;\hat{\psi}_{nk})$, and $\partial_\psi \hat{G}_t(\hat{\theta}_n;\hat{\psi}_{nk})$ with $\hat{G}_t(\hat{\theta}_n;\hat{\psi}_{nk}) = \text{vec}[\partial_\theta y_t(\hat{\theta}_n;\hat{\psi}_{nk})]$. Then, compute $\hat{M}_{n} = \frac{1}{n} \sum_{t=1}^n \partial_{\theta} y_t(\hat{\theta}_n;\hat{\psi}_{nk})^\prime W_n \partial_{\theta} y_t(\hat{\theta}_n;\hat{\psi}_{nk}) + \frac{1}{n} \sum_{t=1}^n [\hat{u}_t^\prime W_n] \otimes \partial_\theta \hat{G}_t(\hat{\theta}_n;\hat{\psi}_{nk})$ and $\hat{D}_{n,\theta,\psi}(k) = \frac{1}{n} \sum_{t=1}^n \partial_{\theta} y_t(\hat{\theta}_n;\hat{\psi}_{nk})^\prime W_n \partial_{\psi} y_t(\hat{\theta}_n;\hat{\psi}_{nk}) + \frac{1}{n} \sum_{t=1}^n [\hat{u}_t^\prime W_n] \otimes \partial_\psi \hat{G}_t(\hat{\theta}_n;\hat{\psi}_{nk})$. Using the same Gaussian quasi-Likelihood terms as above, evaluate $\hat{S}_{t,k} = \hat{M}_{n}^{-1} \{ \partial_\theta y_t(\hat{\theta}_n;\hat{\psi}_{nk})^\prime W_n \hat{u}_t + \hat{D}_{n,\theta,\psi}(k) H_n(\hat{\psi}_{nk})^{-1} \partial_{\psi} L_t(\hat{\psi}_{nk})\}$. Finally, $\hat{V}_{n,k}$ is the estimator for the long-run variance of $\hat{S}_{t,k}$. \subsection{Specification Testing} The population loss $Q(\theta_0;\psi_0)$ defines a distance between the VMA($\infty$) representations of $\tilde{y}_t$ and $y_t(\theta;\psi_0)$. When the model is correctly specified in terms of second-order moments, the minimizer $\theta_0$ yields $Q(\theta_0;\psi_0)= 0$. When the model is misspecified, however, the minimum is strictly positive: $Q(\theta_0;\psi_0) > 0$. The following considers a specification test based on the sample analog $Q_n(\hat{\theta}_n;\hat{\psi}_{nk})$ of the optimal transport distance $Q(\theta_0;\psi_0)$. \begin{assumption} Suppose that: (i). $[ k \log(n)]^8/n = o(1)$; (ii). $[\log(n)]^{4}/k = o(1)$; (iii). $\mathbb{E}(\|e_t\|^{16}) < \infty$; (iv). $\alpha(j) \leq C (1+j)^{-(a+\varepsilon)}$ for $a \geq 6$, $\varepsilon > 0$ and all $j \geq 1$; (v). Assumption (ref) holds with $b \geq 6$; (vi). $\|W_n - W\| = O_p( n^{-1/2} )$, (vii) $\sqrt{nk} \sum_{j=k+1}^\infty \|\Psi_j\|_{op} = o([\log(n)]^{-2})$. \end{assumption} Assumption (ref) is more restrictive than those needed for Theorems (ref) and (ref). When the model is correctly specified, the loss $Q_n(\hat{\theta}_n;\hat{\psi}_{nk})$ is asymptotically determined by the distance between $\hat{\psi}_{nk}$ and $\psi_k$. The distributional results below build on a strong approximation result for NED processes with dependence changing with the lag structure, indexed by $k$. \begin{theorem}[Specification Test] Suppose the conditions for Theorems (ref) and (ref) and Assumption (ref) hold. If the model is correctly specified, then: \begin{align*} nQ_n(\hat{\theta}_n;\hat{\psi}_{nk}) &= n \mathcal{Z}_{n,k}^\prime M_k \mathcal{Z}_{n,k} + o_p(k^{1/2}[\log(n)]^{-2}), \end{align*} where $M_k = \mathbb{E} \left[ (\partial_\psi y_t(\theta_0;\psi_{k0}) + \partial_\theta y_t(\theta_0;\psi_k) M^{-1}E_k )^\prime W (\partial_\psi y_t(\theta_0;\psi_{k0}) + \partial_\theta y_t(\theta_0;\psi_k) M^{-1}E_k ) \right]$ , with $E_k = - \mathbb{E}[ \partial_\theta y_t(\theta_0;\psi_{k0})^\prime W \partial_\psi y_t(\theta_0;\psi_{k0}) ]$, $\mathcal{Z}_{n,k} \sim \mathcal{N}( 0, S_{n,k}/n )$, $S_{n,k} = n\text{var}[\overline{Z}_{n,k}]$, and $M$, $\overline{Z}_{n,k}$ are defined in Theorem (ref). If $S_{n,k}$ and $M_k$ are such that $\text{trace}\left( S_{n,k}M_k \right) \geq O(k)$ and $\text{trace}\left( [S_{n,k}M_k]^2 \right) \geq O(k)$, then for any $\alpha \in (0,1)$: \[ \mathbb{P}\left( nQ_n(\hat{\theta}_n;\hat{\psi}_{nk}) > c_{n,k}(1-\alpha) \right) = \alpha + o(1), \] where $c_{n,k}(1-\alpha)$ is the $1-\alpha$ quantile of $n \mathcal{Z}_{n,k}^\prime M_k \mathcal{Z}_{n,k}$. \end{theorem} Theorem (ref) shows that $nQ_n(\hat{\theta}_n;\hat{\psi}_{nk})$ can be approximated by a weighted sum of independent $\chi^2_1$ random variables. The derivatives $\partial_\psi y_t(\theta;\psi_{k0})$ are given in Lemma (ref). lutkepohl2005 provides formulas for $S_{n,k}$ in the homoskedastic case. The conditions $\text{trace}(S_{n,k}M_k) \geq O(k)$ and $\text{trace}([S_{n,k}M_k]^2) \geq O(k)$ are analogous to the rank conditions needed to ensure the J-test for GMM has a $\chi^2_{k-d}$ distribution, where $k$ is the number of moments and $d$ the number of parameters. The Theorem states that under the null hypothesis of correct specification, the asymptotic size of the test is $\alpha$. The test is also consistent against distant alternatives, see Lemma (ref) in Appendix (ref). Power against local alternatives depends on the ratio $k/n$. A detailed analysis of local power is left to future research. The test in Theorem (ref) involves all variables $\tilde{y}_t$ used in the estimation. In certain settings, the researcher might inquire how well the model fits a specific variable $\tilde{y}_{t,j}$, e.g. consumption if the object of interest is welfare. The following Corollary specializes to a single variable, using a selection matrix $D_j$. The proof is the same as Theorem (ref), and it is omitted. \begin{corollary}[Specification Test on a Single Variable] Suppose the conditions for Theorems (ref) and (ref) and Assumption (ref) hold. Let $Q_{n,j}(\hat{\theta}_n;\hat{\psi}_{nk}) = \frac{1}{n} \sum_{t=1}^n \| y_{t} - \tilde{y}_{t}\|^2_{D_j W_n D_j}$ for $j \in \{1,\dots,d\}$, where $D_j = \text{diag}(\mathbbm{1}_{j=1},\dots,\mathbbm{1}_{j=d})$. If the model is correctly specified, then: \[ nQ_{n,j}(\hat{\theta}_n;\hat{\psi}_{nk}) = n \mathcal{Z}_{n,k}^\prime M_{k,j} \mathcal{Z}_{n,k} + o_p(k^{1/2}[\log(n)]^{-2}),\] where $M_{k,j} = \mathbb{E} \left[ (\partial_\psi y_t(\theta_0;\psi_{k0}) + \partial_\theta y_t(\theta_0;\psi_k) M^{-1}E_{k,j} )^\prime D_j W D_j (\partial_\psi y_t(\theta_0;\psi_{k0}) + \partial_\theta y_t(\theta_0;\psi_k) M^{-1}E_k ) \right],$ with $E_{k,j} = - \mathbb{E}[ \partial_\theta y_t(\theta_0;\psi_{k0})^\prime D_j W D_j \partial_\psi y_t(\theta_0;\psi_{k0}) ]$, and $M$ and $\mathcal{Z}_{n,k}$ are as defined in Theorem (ref). If $S_{n,k}$ and $M_{k,j}$ are such that $\text{trace}\left( S_{n,k}M_{k,j} \right) \geq O(k)$ and $\text{trace}\left( [S_{n,k}M_{k,j}]^2 \right) \geq O(k)$, then for any $\alpha \in (0,1)$: $\mathbb{P}\left( nQ_{n,j}(\hat{\theta}_n;\hat{\psi}_{nk}) > c_{n,k,j}(1-\alpha) \right) = \alpha + o(1)$, where $c_{n,k,j}(1-\alpha)$ is the $1-\alpha$ quantile of $n \mathcal{Z}_{n,k}^\prime M_{k,j} \mathcal{Z}_{n,k}$. \end{corollary} \section{Monte Carlo Simulations} The Monte Carlo simulations are based on the lubik2004 model, given by ((ref)). The observables are log levels of output, inflation, and interest rate (both annualized), which satisfy $ Y_{t}=(0,$ $\pi ^{\ast },$ $\pi ^{\ast }+r^{\ast })^{\prime }+(y_{t},$ $4\pi _{t},$ $4r_{t})^{\prime }$, where output is detrended, and $\pi ^{\ast }$ and $r^{\ast }$ are annualized steady-state rates of inflation and real interest rate with $\beta =(1+r^{\ast }/100)^{-1/4}$. The data are generated using the posterior means from Bayesian inference on the full sample with LS's prior.\footnote{Table (ref), Appendix (ref), includes a description of the parameters, the bounds imposed on the parameters, and the prior $\pi$ used to regulate the estimates.} The VAR includes a constant and $4$ lags as regressors. \begin{table}[ht] \caption{ LS Model: Average Estimate, Standard Deviation, Rejection Rates } \setlength\tabcolsep{4.0pt} { \begin{tabular}{l|ccccccccccccc} \hline \hline & ${\tau }^{-1}$ & ${r}^{\ast }$ & ${\kappa }$ & ${\psi }_{1}$ & ${\psi }_{2}$ & ${\rho }_{r}$ & ${\rho }_{g}$ & ${\rho }_{z}$ & ${\sigma }_{r}$ & ${\sigma }_{g}$ & ${\sigma }_{z}$ & ${\rho }_{gz}$ & ${\pi }^{\ast }$ \\ \hline \textsc{true} & 3.18 & 1.87 & 0.50 & 1.33 & 0.21 & 0.76 & 0.89 & 0.86 & 0.26 & 0.13 & 0.97 & 0.80 & 4.01\\ \hline & \multicolumn{13}{c}{$n = 192$} \\ \hline \textsc{mean} & 1.93 & 1.84 & 0.39 & 1.29 & 0.18 & 0.73 & 0.86 & 0.82 & 0.26 & 0.20 & 1.04 & 0.60 & 3.94\\ \textsc{std} & 0.42 & 0.33 & 0.13 & 0.16 & 0.02 & 0.05 & 0.04 & 0.05 & 0.04 & 0.03 & 0.17 & 0.21 & 0.63 \\ \hline \textsc{rej}$_c$ & 0.14 & 0.07 & 0.05 & 0.00 & 0.00 & 0.01 & 0.03 & 0.04 & 0.01 & 0.05 & 0.03 & 0.04 & 0.11 \\ \textsc{rej}$_r$ & 0.01 & 0.07 & 0.05 & 0.09 & 0.00 & 0.00 & 0.07 & 0.08 & 0.00 & 0.07 & 0.03 & 0.02 & 0.10 \\ \hline \textsc{len}$_c$ & 5.16 & 1.20 & 0.99 & 2.51 & 4.14 & 0.30 & 0.16 & 0.19 & 0.20 & 0.22 & 0.68 & 1.22 & 2.13 \\ \textsc{len}$_r$ & 7.66 & 1.25 & 1.99 & 6.29 & 11.21 & 0.78 & 0.18 & 0.24 & 0.35 & 0.40 & 0.95 & 2.25 & 2.17 \\ \hline \hline & \multicolumn{13}{c}{$n = 500$} \\ \hline \hline \textsc{mean} & 2.08 & 1.87 & 0.38 & 1.30 & 0.18 & 0.75 & 0.89 & 0.84 & 0.26 & 0.16 & 1.03 & 0.63 & 4.08\\ \textsc{std} & 0.45 & 0.23 & 0.11 & 0.11 & 0.02 & 0.03 & 0.02 & 0.03 & 0.03 & 0.02 & 0.11 & 0.15 & 0.44 \\ \hline \textsc{rej}$_c$ & 0.25 & 0.07 & 0.14 & 0.00 & 0.00 & 0.00 & 0.01 & 0.07 & 0.04 & 0.04 & 0.04 & 0.01 & 0.10 \\ \textsc{rej}$_r$ & 0.04 & 0.07 & 0.03 & 0.00 & 0.00 & 0.00 & 0.02 & 0.06 & 0.01 & 0.06 & 0.04 & 0.01 & 0.10 \\ \hline \textsc{len}$_c$ & 3.59 & 0.81 & 0.65 & 1.44 & 2.67 & 0.18 & 0.10 & 0.10 & 0.14 & 0.12 & 0.43 & 0.84 & 1.57\\ \textsc{len}$_r$ & 4.50 & 0.82 & 0.91 & 2.84 & 5.96 & 0.37 & 0.10 & 0.12 & 0.18 & 0.16 & 0.50 & 1.17 & 1.58\\ \hline \hline \end{tabular} \linespread{0.2} \captionsetup{justification=justified} \caption*{ \textbf{Legend:} 200 Monte Carlo replications. \textsc{mean}/\textsc{std}: average and empirical standard error of estimates. \textsc{rej}$_c$, \textsc{rej}$_r$: rejection rates for $5\%$ level t-test. \textsc{len}: median length of $95\%$ confidence intervals. } } \end{table} The baseline sample size corresponds to the full sample estimation below with $n=192$. A larger sample size of $n = 500$ is also considered. The prior from LS is used to regularize the estimates. Table (ref) reports the averages and standard deviations of the estimates, rejection rates using standard errors that assume correct specification and those that allow for misspecification, and the length of resulting $95\%$ confidence intervals. The standard error estimates used to compute the tests and confidence intervals do not account for the prior regularization, which is assumed to be \textit{asymptotically} negligible. Most estimates are centered at the true value when $n=192$. A few estimates are somewhat biased towards the prior mode, most notably the risk aversion $\tau^{-1}$, which lies between the true value $3.18$ and the prior mode $1.88$. The rejection rates are generally close to the 5% level or conservative; significant overrejection is observed only for $\tau^{-1}$ when using the non-robust standard errors, driven by prior-induced bias. The robust standard errors tend to be larger, producing lower rejection rates and wider confidence intervals. The estimation precision improves when $n$ is increased to $500$, and the other conclusions remain similar. The specification test for all variables has a rejection rate of 0.05 and 0.04 for $n=192$ and $n = 500$, respectively. For consumption only, the rejection rates are $0.05$ and $0.04$. Both are close to the nominal level. Additional Monte Carlo simulation for the medium-scale smets2007 model can be found in Table (ref), Appendix (ref). \section{Empirical Illustrations} Three macroeconomic and one financial applications illustrate different aspects of the OT filter and estimation. The first two revisit a small and medium-scale DSGE model using the same sample period 1960Q1-2007Q4 for both. \subsection{Small New-Keynesian Model} This empirical application further considers the LS model. Following LS, there are two specifications: determinacy, with a unique equilibrium, and indeterminacy, where sunspot equilibria exist. The parameters are described in Table (ref). The sample is constructed and divided into subsamples as in clarida2000: the full sample (1960Q1-2007Q4), the pre-Volcker period (1960Q1-1979Q2), and the post-Volcker period (1979Q3-2007Q4). They are associated with determinate, indeterminate, and determinate policy regimes, respectively. To remain consistent with LS and subsequent analyses, their log-prior density $\pi$ is used to penalize the OT loss so that the estimates minimize $n Q_n(\theta;\hat{\psi_n}) - \log \pi(\theta)$. The baseline auxiliary model is a VAR($4$), results with $k=2$ are reported in Appendix (ref). The weighting matrix $W_n$ is diagonal with the inverse of the variances of the three observables. Point estimates and standard errors are reported in Table (ref). The estimates for the pre-Volcker and post-Volcker periods are in line with those in LS, computed using Bayesian likelihood inference. In constrast to likelihood-based estimation, the methods developped here enable us to contrast the actual data with their model-implied values (i.e. the coupling) to obtain an intuitive understanding of their discrepancies. Figure (ref) contrasts actual and model-consistent data for GDP, inflation, and interest rate series, respectively, for the full sample and the two subsamples. For GDP, the actual data exhibit a deeper recession and lower inflation rates in the 1980s than those implied by the model. In other words, the model overpredicts the levels of GDP and inflation compared to the data. This finding confirms that this model, with time-invariant parameters, is unable to capture the rich GDP and inflation dynamics present in the data for the full sample period. \begin{figure}[ht] \caption{LS Model: Actual and Fitted Values} \linespread{0.2} \captionsetup{justification=justified} \caption*{ \textbf{Legend:} Black solid line = data, Purple solid line = coupling. $R^2$ is computed for each variable. } \end{figure} The formal specification test, in Table (ref), investigates tensions between data and model. In the full sample, the test rejects inflation at the 5% significance level. With $k=2$ the model is rejected overall and for each variable individually (Table (ref), Appendix (ref)). The pre and post-Volcker samples do not reject the model. This finding is consistent with LS. None of the variables are individually rejected on the two subsamples. \begin{table}[ht] \caption{LS Model: Specification Test ($k = 4$ lags)} \setlength\tabcolsep{4.0pt} { \ \begin{tabular}{ll|ccc|ccc|ccc} \hline\hline \multicolumn{2}{c|}{\multirow{2}{*}{Specification Test}} & \multicolumn{3}{c|}{Determinacy} & \multicolumn{3}{c|}{Indeterminacy} & \multicolumn{3}{c}{Determinacy} \\ & & \multicolumn{3}{c|}{(Full Sample)} & \multicolumn{3}{c|}{(Pre-Volcker)} & \multicolumn{3}{c}{(Post-Volcker)} \\ \hline \multicolumn{2}{c|} & \textsc{stat} & $10\%$ & $5\%$ & \textsc{stat} & $10\%$ & $5\%$ & \textsc{stat} & $10\%$ & $5\%$ \\ \hline \multicolumn{2}{c|}{All} & 121.2 & 140.0 & 189.6 & 59.4 & 106.2 & 169.6 & 65.9 & 164.8 & 235.3 \\ \hline \multicolumn{2}{c|}{Output} & 56.0 & 90.6 & 123.8 & 42.8 & 51.1 & 86.0 & 30.1 & 131.8 & 188.6 \\ \multicolumn{2}{c|}{Inflation} & 45.2 & 27.1 & 37.8 & 7.0 & 33.9 & 55.0 & 16.6 & 22.2 & 29.0 \\ \multicolumn{2}{c|}{Interest Rate} & 20.0 & 31.1 & 42.1 & 9.6 & 22.7 & 36.6 & 19.2 & 22.0 & 30.8 \\ \hline\hline \end{tabular} \linespread{0.2} \captionsetup{justification=justified} \caption*{ \textbf{Legend:} \textsc{stat}: $nQ_n(\hat{\theta}_n;\hat{\psi}_{nk})$. $10\%$, $5\%$: critical values for specification test at corresponding significance levels. All: joint test on all variables. Output, Inflation, Interest Rate: test on individual variables. $n = 192,78,114$ for the full, pre and post-Volcker samples. } } \end{table} \subsection{Medium-Scale DSGE Model} The second empirical application considers the smets2007 model. Table (ref), Appendix (ref), includes parameter interpretations and prior distributions from SW used in our estimation. This model includes 36 free parameters and is estimated on 7 observables: consumption, investment, output and wage growth, hours worked, inflation, and interest rate. There are as many shocks: productivity, exogenous spending, monetary policy, investment-specific technology, price markup, wage markup, and risk premium shocks. \begin{table}[htb!] \caption{SW Model: Estimates and Standard Errors (selected coefficients)} \setlength\tabcolsep{4.5pt} { \begin{tabular}{ll|ccc||cc} \hline\hline \multirow{2}{*}{${\theta }$} & \multirow{2}{*}{Parameter Interpretation} & \multicolumn{3}{c||}{OT Estimate} & \multicolumn{2}{c}{Posterior} \\ & & {\ \textsc{est}} & {\ \textsc{sd}$_{c}$} & {\ \textsc{sd}$_{r}$} & {\ \textsc{mean}} & \textsc{std} \\ \hline ${\varphi }$ & {\ Investment adjustment cost} & 3.00 & 1.50 & 3.42 & {\ 6.12} & 0.99 \\ ${\sigma }_{c}$ & {\ Elast. of Intertemporal substitution} & 1.01 & 0.12 & 0.52 & {\ 1.50} & 0.14 \\ ${\lambda }$ & {\ Habit persistence} & 0.74 & 0.12 & 0.38 & {\ 0.71} & 0.04 \\ ${\sigma }_{l}$ & {\ Labor supply elasticity} & 1.00 & 2.38 & 2.92 & {\ 2.25} & 0.55 \\ ${r}_{\pi }$ & {\ Taylor rule: inflation weight} & 1.80 & 1.06 & 0.77 & {\ 2.03} & 0.16\\ ${r}_{\Delta y}$ & {\ Taylor rule: output gap change weight} & 0.16 & 0.05 & 0.10 & {\ 0.21} & 0.02 \\ ${r}_{y}$ & {\ Taylor rule: output gap weight} & 0.16 & 0.16 & 0.08 & {\ 0.10 } & 0.02 \\ ${\rho }$ & {\ Taylor rule: interest rate smoothing} & 0.90 & 0.08 & 0.07 & { \ 0.82} & 0.02 \\ ${\rho }_{a}$ & {\ Productivity shock AR} & 0.94 & 0.03 & 0.22 & {\ 0.97} & 0.01 \\ ${\rho }_{b}$ & {\ Risk premium shock AR} & 0.68 & 0.10 & 0.19 & {\ 0.28} & 0.12 \\ ${\rho }_{g}$ & {\ Exogenous spending shock AR } & 0.86 & 0.09 & 1.22 & {\ 0.97} & 0.01 \\ ${\rho }_{i}$ & {\ Investment shock AR} & 0.47 & 0.15 & 0.13 & {\ 0.70} & 0.06 \\ ${\rho }_{r}$ & {\ Monetary policy shock AR} & 0.47 & 0.28 & 0.71 & {\ 0.17} & 0.07 \\ \hline\hline \end{tabular} \linespread{0.2} \captionsetup{justification=justified} \caption*{ \textbf{Legend:} Full table of estimates, prior distribution, and estimation bounds can be found in Tables (ref), (ref), Appendix (ref). \textsc{mean}, \textsc{std} = posterior mean and standard deviation. } } \end{table} \textbf{Estimation.} To illustrate the scope of OTE, the full model and $3$ singular versions are estimated using the same method. Estimates are reported in Tables (ref) (for the full model) and (ref) (for singular models), Appendix (ref). Likelihood-based posterior estimates computed with SW's prior are reported in Table (ref) as a reference. All specifications rely on a VAR(4) auxiliary model; the inverse of the variances of the observables is used as $W_n$. The $3$ singular models remove, in order, the risk premium, wage markup, and price markup shocks reducing to $6$, $5$, and $4$ shocks for $7$ observables. The choice of shocks to remove follows qu2018, to reflect a view that they have a weaker structural interpretation than the remaining 4. It is interesting to examine their impact on the model fit. The standard likelihood approach cannot estimate singular DSGE models because the covariance matrix of the one-step-ahead forecasting errors $\Sigma(\theta)$ is singular. qu2018 used a composite likelihood framework and did not formally test the resulting models. OTE handles both singular and nonsingular models within the same framework. Additionally, note that although SW chose to fit their model to seven variables, it has implications for additional macro variables including the price of capital and capital utilization rate. Adding any of these variables to the set of observables will immediately make the model singular. In fact, for all medium-scale DSGE models, nonsingularity arises only because we restrict the estimation to a limited set of macro variables. For the full model, the OT estimates are similar to the posterior means. In all but two cases, the posterior means fall within the 95% confidence intervals obtained from OT estimates and robust standard errors. For the remaining two cases, OT produces a more persistent risk premium shock process with a lower residual standard deviation. The robust standard errors are almost always greater than those assuming correct model specification. For singular models, removing the risk premium shock has little effect on the estimates --- they are close to the nonsingular case; all confidence intervals overlap with their nonsingular counterparts. When the wage markup shock is also removed, ${\small \xi }_{w}$ (wage stickiness) and $\iota _{w}$ (wage indexation) decrease, while $\rho _{r}$ (monetary policy shock persistence) increases, though their confidence interval still overlaps due to substantial estimation uncertainty. When further removing the price markup shock, ${\small \xi }_{p}$ (price stickiness) and ${\small \iota }_{p}$ (wage stickiness) both drop noticeably. Although not reported here, the effects of these parameter differences on the model can be further assessed by plotting the impulse response functions. \paragraph{Specification Testing.} The specification test does not reject the original model at the 5% significance level, as shown in Table (ref). However, the test rejects the model's fit for consumption, even with 7 shocks. Further investigation reveals that the model under-predicts contractions, e.g. in 1974Q4 consumption fell by 2.38% vs. 1.34% for the fitted values. Fitted consumption is less volatile (standard deviation of 0.53 vs. 0.68), and more persistent (autocorrelation of 0.46 vs. 0.18). For singular models, when the risk premium shock is removed, the specification test rejects the full model at the 5% significance level; for individual tests, the results for consumption, wage, and interest rate reject the null hypothesis. When the wage markup shock is removed, the tests on output and labor also reject, implying that 5 out of 7 variables are now rejected. Finally, when the price markup shock is removed, the conclusions remain the same as in the five-shock case, with only two variables—investment and inflation—not rejected by the test at the 5% level. This is the first attempt to formally test singular DSGE models. The results pinpoint model features that remain compatible or become incompatible with data once shocks are removed. They can be useful tools for researchers to determine which latent processes and mechanisms contribute to the fit of a model within a unified framework. \begin{table}[ht] \caption{SW Model: Specification Testing With(out) Stochastic Singularity} \setlength\tabcolsep{4.0pt} { \begin{tabular}{l|ccc||ccc||ccc||ccc} \hline\hline & \multicolumn{3}{c||}{7 shocks} & \multicolumn{3}{c||}{6 shocks} & \multicolumn{3}{c||}{5 shocks} & \multicolumn{3}{c}{4 shocks} \\ & \textsc{stat} & 10% & 5% & \textsc{stat} & 10% & 5% & \textsc{stat} & 10% & 5% & \textsc{stat} & 10% & 5% \\ \hline All & 256.6 & 320.7 & 395.0 & 411.4 & 282.1 & 364.5 & 597.7 & 291.4 & 357.0 & 632.0 & 288.6 & 353.3 \\ \hline Cons. & 60.4 & 31.9 & 36.8 & 56.5 & 32.6 & 39.9 & 192.9 & 6.3 & 8.9 & 192.7 & 6.1 & 8.5 \\ Invest. & 25.6 & 33.2 & 38.6 & 34.1 & 31.9 & 36.6 & 45.6 & 41.0 & 47.3 & 45.0 & 41.1 & 47.7 \\ Output & 34.3 & 37.0 & 41.8 & 34.9 & 32.4 & 36.7 & 54.4 & 39.8 & 45.1 & 57.9 & 39.2 & 44.8 \\ Labor & 16.3 & 35.4 & 47.2 & 14.7 & 48.4 & 66.8 & 86.9 & 49.3 & 63.5 & 97.6 & 47.8 & 60.0 \\ Infl. & 52.8 & 92.5 & 122.7 & 75.1 & 93.5 & 124.2 & 67.6 & 101.1 & 137.8 & 69.6 & 104.1 & 146.1 \\ Wage & 22.6 & 64.8 & 80.3 & 108.3 & 26.2 & 33.3 & 76.4 & 30.9 & 37.6 & 98.3 & 27.4 & 33.1 \\ Int. Rate & 44.6 & 47.9 & 63.3 & 87.9 & 33.8 & 45.2 & 73.9 & 49.9 & 67.3 & 70.9 & 48.8 & 65.6 \\ \hline\hline \end{tabular} \linespread{0.2} \captionsetup{justification=justified} \caption*{ \textbf{Legend:} All: specification test on all $7$ variables (consumption, investment, output, labor, inflation, wage, interest rate). \textsc{stat}: test statistic for specification test. $5\%$, $10\%$: critical values.} } \end{table} \textbf{Filtering the Latent Shock Processes.} The OTF produces model-consistent values of latent variables, including shock processes. Using the original SW model, we compare these filtered values with their counterparts produced by the KF, which does not enforce model consistency. The same parameter values (OT estimates) are used to ensure comparability. \begin{figure}[h] \caption{SW Model: Filtered Shock Processes} \\ \linespread{0.2} \captionsetup{justification=justified} \caption*{ \textbf{Legend:} Black solid line: Kalman Filter (KF), Purple solid line: Optimal Transport Filter (OTF). Both filters are applied using the same OT estimates found in Table (ref), Appendix (ref). Orange line: Real Government Consumption Expenditures and Gross Investment; Source: FRED (GCEC1). } \end{figure} Figure (ref), panel a) displays the results for the 7 shocks separately. For TFP, the KF yields a puzzling conclusion: the economy was under a positive TFP process from 1960 until about 1980, and then affected by a negative TFP process between 1980 and 2000. In contrast, the OT filter produces a process with negative values during the early 1960s, the mid-1970s recessions, and the slowdown leading up to the 2001 recession. These estimates are clearly more interpretable than the KF values. The investment technology shock shows a similar pattern: the KF yields a mostly positive process from the 1970s to the mid-1980s, which then switches to a mostly negative process. The OT estimates do not have this problem. For exogeneous spending, KF yields negative values most of the time and exhibits a downward trend, in sharp contrast with the zero-mean assumption. The OTF estimates do not have this problem. For the monetary policy shock, KF produces large oscillating values in the early 1980s, which is puzzling given the monetary tightening that characterizes this period. In contrast, the OTF produces mostly positive values for this period, consistent with this characterization. The last $3$ shocks show closer resemblances between KF and OTF. Table (ref), Appendix (ref), presents the cross and serial-correlations of the 7 shock processes. OT consistently gives values close to the true ones, while KF shows significant discrepancies in several cases; for example, the true first-order serial correlations for investment and monetary policy shocks are both 0.47, but KF yields 0.83 and 0.09. This illustrates that, when using KF, the parameter estimates may not capture the dynamics of the filtered values. As a validation exercise, Figure (ref), panel b), compares the exogenous spending shocks process with government consumption data. They are closely related since exogeneous spending is the difference between output and the sum of consumption, investment and capital utilization smets2007. The correlation between the data and filtered values is 0.36 and 0.05 for OTF and KF, respectively. Although this data was \textit{not used} in the estimation, the OTF captures some of its variation using other series and the model. These results demonstrate that the OT filter, by enforcing model consistency, produces filtered values that obey model assumptions and can be more interpretable in practice. To wrap up, we have considered this medium-scale DSGE model to illustrate that the proposed methods can be used to estimate singular and nonsingular models, testing their specifications, and produce filtered variables, all within the same framework. The methods' applications are not restricted to macroeconomics; we next consider a financial application. \subsection{Affine Term Structure Model} The third application considers a term structure model where three latent factors explain six yields. Algorithm (ref) provides a way to assess the extent to which the unaltered, stochastically singular, structural model fits the empirical data; no measurement errors are introduced. The specification of the structural model follows ang2003 and hamilton2012. Let $F_{t}$ denote three latent factors, which follow a Gaussian VAR: \( F_{t+1}=c+\rho F_{t}+v_{t+1}, \) with $v_{t+1} \overset{iid}{\sim} \mathcal{N}(0,I)$. Affine term structure models assume that the one-period short rate $r_{t}$ is an affine function of the state vector. They specify the stochastic discount factor $M_{t,t+1}$ to be a function of $\lambda_{t}$, the market prices of risk, which is also an affine function of $F_{t}$: $r_{t}=\delta _{0}+\delta _{1}^{\prime }F_{t}$, $M_{t,t+1}=\exp [-r_{t}-(1/2)\lambda _{t}^{\prime }\lambda _{t}-\lambda_{t}^{\prime }v_{t+1}]$, $\lambda _{t}=\lambda +\Lambda F_{t}$. As highlighted in hamilton2012, under the risk-neutral measure all assets are discounted by the short-term interest rate: $ P_{t}=\exp (-r_{t})\int P_{t+1}(F_{t+1})\phi (F_{t+1};c^{Q}+\rho ^{Q}F_{t},I)dF_{t+1}$, where $F_{t+1}=c^{Q}+\rho ^{Q}F_{t}+v_{t+1}^{Q}$, with $v_{t+1}^{Q} \overset{iid}{\sim} \mathcal{N}(0,I)$, $c^{Q}=c-\lambda $, and $\rho ^{Q}=\rho -\Lambda $. It is well-documented that parameter normalizations are necessary to identify the model. Following ang2003 and hamilton2012, we set: $c=0,\delta _{1}\geqslant 0$, and $\rho ^{Q}$ lower triangular. Then, the yield $y_{t}^{n}$ on an $n$-period pure-discount bond is given by: \begin{align*} y_{t}^{n} &= a_{n}+b_{n}^{\prime }F_{t}, \quad b_{n} =\frac{1}{n}\left( I+\rho ^{Q\prime }+...+\left( \rho ^{Q\prime }\right) ^{n-1} \right) \delta _{1} \\ a_{n} &=\delta _{0}+\frac{1}{n}\left( b_{1}^{\prime }+2b_{2}^{\prime }+...+(n-1)b_{n-1}^{\prime }\right) c^{Q}-\frac{1}{2n} \left( b_{1}^{\prime }b_{1}+4b_{2}^{\prime }b_{2}+...+(n-1)^{2}b_{n-1}^{\prime }b_{n-1} \right). \end{align*} The parameters to be estimated are $\delta _{0},\delta _{1},\rho ^{Q},c^{Q}$ , and $\rho$. The diagonal elements of $\rho ^{Q}$ are required to be in decreasing order to ensure identification. The auxiliary model is VAR($4$). The yields are weighted diagonally by the inverse of their variance. The OT objective $Q_n$, is penalized by prior distribution which enforces an ordering of the factors. \begin{table}[ht] \caption{Affine Term Structure Model: Parameter Estimates} \setlength\tabcolsep{4.0pt} { \begin{tabular}{c|ccc|ccc} \hline \hline & \multicolumn{3}{c|}{Prior [\textsc{mean}, \textsc{sd}]} & \multicolumn{3}{c}{ OT Estimates} \\ \hline \multirow{3}{*}{$\rho ^{Q}$} & [0.9,0.2] & [0,0.2] & - & 0.999 & - & - \\ & [0,0.2] & [0.8,0.2] & - & 0.022 & 0.963 & - \\ & [0,0.2] & [0,0.2] & [0.6,0.2] & 0.018 & 0.209 & 0.723 \\ \hline $\delta _{1}$ & [0.01,0.1] & [0.01,0.1] & [0.01,0.1] & 0.003 & 0.023 & 0.037 \\ \hline $\delta _{0}$ & [0.4,0.2] & - & - & 0.441 & - & - \\ \hline $c^{Q}$ & [0,1] & [0,1] & [0,1] & 1.254 & -0.026 & 0.536 \\ \hline \multirow{3}{*}{$\rho$} & [0.9,0.2] & [0,0.2] & [0,0.2] & 0.958 & 0.012 & 0.071 \\ & [0,0.2] & [0.8,0.2] & [0,0.2] & 0.004 & 0.908 & 0.102 \\ & [0,0.2] & [0,0.2] & [0.6,0.2] & 0.011 & 0.133 & 0.770 \\ \hline \hline \end{tabular} \linespread{0.2} \captionsetup{justification=justified} \caption*{ \textbf{Legend:} A weakly informative prior is used to enforce the ordering of the factors. } } \end{table} The estimates in Table (ref) are similar to hamilton2012.\footnote{hamilton2012 use four yields and introduce measurement error in one of the series.} The first factor is very persistent, and the off-diagonal elements of $\rho$ are small in magnitude. The estimate of $\delta_{0}$ is close to the mean of the short-term rate. The main difference is that the first element of $\delta _{1}$ in our case is smaller than theirs (their estimate is 0.017). The parameters $\rho^{Q},\delta _{1}$ and $c^{Q}$ are highly correlated: it is possible to move their values jointly with little effect on $Q_n$, suggesting weak identification.For this reason standard errors are not reported. The specification test is not reported either. The prior helps stabilize the estimates. \begin{figure}[ht] \caption{Affine Term Structure Model: Actual and Model-Consistent Yields} \linespread{0.2} \captionsetup{justification=justified} \caption*{ \textbf{Legend:} Actual yields: Black solid line, Coupling: Purple solid line. Fit shown for prior regularized estimates; the fit is virtually identical without prior regularization.} \end{figure} What is intriguing is whether, with three shocks, the model can approximate the dynamics of six observables. Figure (ref) compares actual yield data with their couplings. These values track each other closely, with the $R^{2}$ equal to 0.97 in all six cases. With three shocks, the model accounts for 97% of the dynamics in the data for this sample period. The plots suggest no apparent model misspecification for this sample period. Adding measurement errors in this case would be a shortcut to obtaining parameter estimates with a likelihood and cannot reveal additional information regarding the baseline model specification. Our approach allows us to obtain parameter estimates with a graphical method to assess the model fit. Finally, we conjecture that incorporating more information to improve the identification of $\rho^{Q},\delta _{1}$ and $c^{Q}$ can be beneficial for further improve this model. \subsection{Trend-Cycle Decomposition} The last application illustrates how the methodology can be applied in a setting where the stationarity assumptions used in the theoretical analysis are not satisfied. Trend-cycle decompositions are routinely used to date business cycles for many countries. One approach is to model the log real GDP using an unobserved component model as in watson1986: \[ \underbrace{y_t = \tau_t + c_t}_{\text{ $\log$(GDP) }}, \text{ where } \underbrace{\tau_t = \mu + \tau_{t-1} + \eta_t}_{\text{trend component}} \text{ and } \underbrace{c_t = \rho_1 c_{t-1} + \rho_2 c_{t-2} + e_t}_{\text{cycle component}}, \] where $(\eta_t,e_t) \overset{iid}{\sim} \mathcal{N}(0, \text{diag}(\sigma_{\eta}^2,\sigma_{e}^2) )$. While this DGP is non-stationary, and thus not covered by the parameter estimation results above, it illustrates that the methodology applies to a broader set of filtering problems. \begin{figure}[H] \caption{U.S. GDP -- 1947Q1-2023Q4 -- Trend-Cycle Decomposition} \linespread{0.2} \captionsetup{justification=justified} \caption*{ \textbf{Legend:} Blue dashed line: linear trend w. cosines $\delta_0 + \delta_1 t + \delta_2 \cos( 2 \pi t/n ) + \delta_3 \cos( 2 \pi 2 t/n ) + \delta_4 \cos( 2 \pi 3 t/n ) + \delta_5 \cos( 2 \pi 4 t/n )$. Red dashed line: Kalman Filter estimates of $\tau_t$. Vertical bars: NBER recession dates. } \end{figure} Figure (ref) displays U.S. log-GDP between 1947 and 2023 (left panel) and the extracted cycle components (right panels). The left panel indicates that the KF (red dashed line) does not fully capture changes in the trend growth rate over the three-quarters of a century spanned by the data. In particular, the estimated trend is systematically below the log-GDP between 1965 and 2008. This issue was already raised by perron2009. A flexible trend estimate (blue dashed line), described below the figure, better captures the gradual changes in the trend component. The number of cosine terms is chosen to have the shortest periodicity at 20 years, larger than the business cycle frequency, i.e. 1.5 to 8 years. The top right panel reports the KF estimates for the cycle component. The KF estimated trend is systematically below GDP throughout the period 1965-2008 (left panel). Likewise, the filtered cycle component (top right panel) is systematically positive during this period. This would suggest that the economy is characterized by a continuous expansion over this forty-year period. Clearly, the cycle estimates are contaminated by misspecification in the trend. Finally, the bottom right panel shows the OTF values, computed using the flexible trend as a basis for the reduced form model. See Appendix (ref) for details. No terms were added to handle the Covid recession. Unlike KF estimates, it marks turning points (vertical bars) with good accuracy and visually appears to be stationary with zero mean. \section{Extension: Non-Linear State-Space Models} The following provides a general-purpose Algorithm to perform Optimal Transport Filtering for non-linear non-Gaussian state-space models. It is similar to Algorithm (ref) with the difference that the \textbf{Predict}, \textbf{Transport}, and \textbf{Update} steps accomodate more general state-space models ((ref)). It is not closed form and more computationally demanding than Algorithm (ref). \begin{algorithm}[H] \caption{Optimal Transport Filter} { \begin{algorithmic}[1] \Procedure{\textsc{otf}}\newline \textbf{Inputs:} 1) Sample: $\tilde{y}_1,\dots,\tilde{y}_n$, predictive distribution $\tilde{p}(\tilde{y}_t|\tilde{y}_{t-1},\dots)$\newline\hphantom{\textbf{Inputs:}} 2) Model: conditional distribution $p(y_t,z_t|z_{t-1};\theta)$, initial beliefs $z_0 \sim p_{0|0}(z_0)$\newline \textbf{Outputs:} 1) Mapped data $y_1,\dots,y_n$, 2) Filtered states $z_{t|t} \sim p( z_t|y_t,\dots,y_1 )$. \For{$t \in \{1,\dots,n\}$} \State{\textbf{Predict:} Using the model, compute} \Comment{(Filter)}\newline \hphantom{\textbf{Inputs:}} $p_{t|t-1}(y_t,z_t) = \int p(y_t,z_t|z_{t-1})p_{t-1|t-1}(z_{t-1})dz_{t-1}$ \State{\textbf{Transport:} Find $\pi_{t|t-1}$ which solves $\min_{\pi \in \Pi_{t|t-1}} \mathbb{E}_{\pi}(\|y_t-\tilde{y}_t\|^2)$ } \Comment{(OT)} \newline \hphantom{\textbf{Inputs:}} where $\Pi_{t|t-1} = \{ \pi(y_t,\tilde{y}_t) \text{ s.t. } \int \pi(y_t,\tilde{y}_t)d\tilde{y}_t = p_{t|t-1}(y_t), \int \pi(y_t,\tilde{y}_t)dy_t = \tilde{p}_{t|t-1}(\tilde{y}_t) \}$ \State{\textbf{Update:} $p_{t|t}(z_t) \propto p_{t|t-1}(z_t,y_t)/\pi_{t|t-1}(y_t|\tilde{y}_t)$ } \Comment{(Filter)}\newline \hphantom{\textbf{Inputs:}} using the distributions $p_{t|t-1}$ from step 3 and $\pi_{t|t-1}$ from step 4. \EndFor \EndProcedure \end{algorithmic}} \end{algorithm} The following briefly discusses the main steps of the Algorithm. An estimate of $\tilde{p}_{t|t-1}$ is required as inputs, as for linear state-space models. More details on the filtering (Filter) and optimal transport (OT) steps can be found in Algorithm (ref) in Appendix (ref). \textbf{Predict.} This step can be carried out using simulations as in the particle filter algorithm. For discrete or discretized models, these can be carried out using only matrix operations chopin2020. This produces draws for Transport step below. \textbf{Transport.} When $\tilde{p}_{t|t-1}$ and $p_{t|t-1}$ are non-Gaussian, the transport still has a closed form when both $\tilde{y}_t$ and $y_t$ are scalar: $y_t = F_{t|t-1}^{-1} \circ \tilde{F}_{t|t-1}(\tilde{y}_{t|t-1})$. For multivariate outcomes, the solution is generally not closed form and needs to be computed numerically. From a finite number of draws $\tilde{y}_t^b \sim \tilde{p}_{t|t-1}$ and $y_t^b \sim p_{t|t-1}$, $b=1,\dots,B$, the optimal transport problem can be solved exactly as a linear program or approximated using entropic regularization and Sinkhorn’s algorithm. The latter is implemented in Algorithm (ref). \textbf{Update.} As in Algorithm (ref), the update step is standard but applied to $y_t$ instead of $\tilde{y}_t$ to ensure the filtered variables are model-consistent. \section{Conclusion} This paper has introduced a computationally attractive method for filtering and estimation of potentially misspecified dynamic models using dynamic optimal transportation. Empirical applications illustrate how this can be used to visually assess the fit of a model, by comparing the actual and the coupled time-series, and to formally test the model specification over all or some specific variables. Several extensions could be of interest in future research. Deriving a plugin map for general non-linear state-space models in ((ref)) could be useful to circumvent the curse of dimensionality. A useful Corollary to Theorem (ref) would be to consider an Anderson-Rubin type statistic for models that are potentially weakly identified. \printbibliography[heading=subbibliography]