Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
15,685 characters · 9 sections · 15 citation commands
On Parameter Estimation in Unobserved Components Models subject to Linear Inequality Constraints
Consider the following unobserved components model.
where the initial condition $\tau_0$ is treated as a parameter.
Suppose we wish to restrict the trend $\tau_t$ to be in the interval $[a_t,b_t], t=1,...,T$ where $a_t$ and $b_t$ are constants. Linear inequality constraints like these occur quite naturally in many economic applications. For instance, chan2016bounded discussed such situations in the context of modeling trend inflation. The problem of interest is to efficiently estimate the parameters involved in an unobserved components model under such restrictions. By efficiency, we mean the run time and the sample efficiency.
In the presence of the aforementioned inequality constraints, the conditional distribution of $\boldsymbol{\tau}$ given the other parameters, i.e. $p(\boldsymbol{\tau} | \mathbf{y}, \sigma^2, \omega^2, \tau_0)$ is nonstandard and thus conventional methods of inference in state-space models cannot be used. chan2012estimation suggest an approach for estimation in such non-linear state-space models where the key element is to approximate $p(\boldsymbol{\tau} | \mathbf{y}, \sigma^2, \omega^2, \tau_0)$ using a multivariate Gaussian density. The Gaussian approximation is then used as a proposal density for an acceptance-rejection Metropolis-Hastings algorithm (chan2016bounded).
We propose a new quadratic programming-based approach to approximating the nonstandard density using a Gaussian density. In general, the idea is to efficiently solve the following quadratic programming problem.
The above quadratic program maximizes the kernel of a multivariate Gaussian density for $\boldsymbol{\tau}$ for general linear inequality constraints. In general, we can say the following about the time complexity of solving a quadratic program like the one mentioned above.
As our problem has more structure like $\mathbf{A} = \mathbf{I}_T$. We may be able to further improve on the empirical run time if not the asymptotic polynomial time complexity.
We know the following. $$p(\boldsymbol{\tau} | \mathbf{y}, \sigma^2, \omega^2, \tau_0) \sim \mathcal {N}(\hat{\boldsymbol \tau},\mathbf{K}_{\boldsymbol\tau}^{-1}),$$
where $\mathbf{K}_{\boldsymbol\tau} = \left(\frac{1}{\omega^2} + \frac{1}{\sigma^2} \right)\mathbf{H}'\mathbf{H}$, $\hat{\boldsymbol \tau} = \mathbf{K}_{\boldsymbol\tau}^{-1} \left(\frac{\tau_0}{\omega^2}\mathbf{H}'\mathbf{H} + \frac{1}{\sigma^2}\mathbf{H}'\mathbf{H}\mathbf{y}\right)$, and $\mathbf{H}$ is given as follows.
We set $\mathbf{m}=\hat{\boldsymbol \tau}$, $\mathbf{C}=\mathbf{K}_{\boldsymbol\tau}$, and $\mathbf{A} = \mathbf{I}_T$. We use quadprogram available in MATLAB to solve the following quadratic programming problem.
We can write it equivalently as follows.
We use the US quarterly CPI data from 1947Q1 to 2015Q4. Specifically, given the quarterly CPI figures $z_t$, we compute $y_t = 400(\log z_t - \log z_{t-1})$ and use it as the inflation rate. We take subsets of $\mathbf{m}$ and $\mathbf{C}$ to investigate the run times for different values of $T$.
The run times in solving the above quadratic programming problem for different values of $T$ are reported in Table 1.
Comment. We observed that the empirical run times are quite small. This suggests that the idea is promising and we should explore it further. The solution to the discussed quadratic program is bound to lie within the constraints. Hence, for Bayesian computations, it is a better candidate for the mean of the distribution which the parameter is getting generated from. The idea is to improve the estimation and avoid any approximations due to the linear constraints and also expect that the sample efficiency will improve. The estimation performed using this idea will be termed QuadProg henceforth.
chan2016bounded discuss a new model of trend inflation with inequality constraints. The model is discussed in great detail is the original paper. Let $\pi_t$ is an observed measure of inflation at time $t$, then the measurement and state equations of the model are given as follows.
where $\epsilon_t \sim \mathcal{N}(0,1)$ and $\epsilon_t^h \sim \mathcal{N}(0,\sigma^2)$. The linear inequality constraints on the behavior of $\tau_t$ an $\rho_t$ are imposed by considering the following.
where $\mathcal{TN}(a,b;\mu,\sigma^2)$ denotes the Gaussian distribution with mean $\mu$ and variance $\sigma^2$ truncated to the interval $(a, b)$.
Furthermore, the model considers the following.
where $a_\tau,b_\tau,\tau_0,\omega_\tau^2,\rho_0,\omega_\rho^2,h_0$ and $\omega_h^2$ are known constants.
Following are the priors for the unknown parameters.
chan2016bounded discusses given an initial value $y_0$ of the CPI inflation rate, use the measurement equation to generate the CPI inflation rate data, $\boldsymbol{y}$, and develop a Markov chain Monte Carlo (MCMC) algorithm for generating from the posterior densities for the above model. Due to the presence of the inequality constraints, $p(\boldsymbol{\tau} | \mathbf{y}, \boldsymbol{\rho}, \boldsymbol{h},\boldsymbol{\theta})$ and $p(\boldsymbol{\rho} | \mathbf{y}, \boldsymbol{\tau}, \boldsymbol{h},\boldsymbol{\theta})$ are nonstandard and thus conventional methods of inference in state-space models cannot be used. Instead, they use an approach developed in chan2012estimation for posterior sampling in nonlinear state-space models. An essential element of this algorithm is a Gaussian approximation to the conditional density $p(\boldsymbol{\tau} | \mathbf{y}, \boldsymbol{\rho}, \boldsymbol{h},\boldsymbol{\theta})$. This uses a precision-based algorithm that builds upon results derived for the linear Gaussian state-space model by chan2009efficient. The Gaussian approximation is then used as a proposal density for an ARMH step. They also use this algorithm to draw from $p(\boldsymbol{\rho} | \mathbf{y}, \boldsymbol{\tau}, \boldsymbol{h},\boldsymbol{\theta})$ and $p(\boldsymbol{h} | \mathbf{y}, \boldsymbol{\tau}, \boldsymbol{\rho},\boldsymbol{\theta})$.
We use the quadratic-programming-based approximation of the Gaussian density (QuadProg) as proposed and discussed in this report for modifying the posterior sampling procedure in chan2016bounded for $\boldsymbol{\tau}$ and $\boldsymbol{\rho}$ as follows.
We used the US quarterly CPI data from 1947Q1 to 2011Q3. Specifically, given the quarterly CPI figures $z_t$, we computed $y_t = 400(\log z_t - \log z_{t-1})$ and used it as the inflation rate. As suggested and discussed in chan2016bounded, we set the hyperparameters as follows.
We used 10,000 Monte-Carlo simulations with 1000 as the burn-in time.
We used the following two estimation methods.
We used Matlab for the purpose of implementation.
We estimated the trend for the US quarterly CPI data using AR-Trend-Bound and Quad-Prog methods.
Figure (ref) shows the observed rate of inflation, trend estimated using AR-Trend-Bound, and trend estimated using Quad-Prog. We observe that the trend estimates using the two different methods are quite close to each other confirming that the the proposed modification indeed does good estimation.
One of the limitations of Monte Carlo methods is that different draws are seemingly/pseudo-random, i.e. there may exist a significant autocorrelation in the draws generated over different time points. We are interested in understanding how does the proposed method compares to the existing approach in terms of presence of such autocorrelation. A common diagnostic of sample efficiency is the inefficiency factor, defined as follows.
where $\phi_l$ is the sample autocorrelation at lag length $l$ and $L$ is chosen large enough so that the autocorrelation tapers off. Note that purely random/independent draws would give an inefficiency factor of 1. Inefficiency factors indicate how many extra draws need to be taken to give results equivalent to independent draws. For instance, if we take 5000 draws of a parameter and find an inefficiency factor of 10, then these draws are equivalent to 500 independent draws from a given distribution.
For any given parameter vector of length $T$, there are $T$ inefficiency factors (one for each $t$). As $\boldsymbol{\tau}, \boldsymbol{\rho}$ and $\boldsymbol{h}$ each is of length $T$, we have $3T$ total number of inefficiency factors corresponding each approach of estimation. We used $L=100$ to calculate the inefficiency factors. For reporting, we construct a box plot of $T$ inefficiency factors for each parameter. We compare the box plots for each parameter across different methods of estimation.
Figure (ref) shows the box plots of inefficiency factors for sampling $\boldsymbol{\tau}$ using AR-Trend-Bound and Quad-Prog methods. We observe that Quad-Prog leads to a lower average and maximum sample inefficiency as compared to AR-Trend-Bound.
Figure (ref) shows the box plots of inefficiency factors for sampling $\boldsymbol{\rho}$ using AR-Trend-Bound and Quad-Prog methods. We observe that Quad-Prog leads to a lower average and maximum sample inefficiency as compared to AR-Trend-Bound.
Figure (ref) shows the box plots of inefficiency factors for sampling $\boldsymbol{h}$ using AR-Trend-Bound and Quad-Prog methods. We observe that Quad-Prog AR-Trend-Bound perform very similar in terms of average and maximum sample efficiency.